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Abstract. In two preceding papers HI |2l we have shown that, when reaction networks are 
well-removed from equilibrium, explicit asymptotic and quasi-steady-state approximations 
can give algebraically-stabilized integration schemes that rival standard implicit methods 
in accuracy and speed for extremely stiff systems. However, we also showed that these 
explicit methods remain accurate but are no longer competitive in speed as the network 
approaches equilibrium. In this paper we analyze this failure and show that it is associated 
with the presence of fast equilibration timescales that neither asymptotic nor quasi-steady- 
state approximations are able to remove efficiently from the numerical integration. Based 
on this understanding, we develop a partial equilibrium method to deal effectively with the 
approach to equilibrium and show that explicit asymptotic methods, combined with the new 
partial equilibrium methods, give an integration scheme that plausibly can deal with the stiffest 
networks, even in the approach to equilibrium, with accuracy and speed competitive with that 
of implicit methods. Thus we demonstrate that such explicit methods may offer alternatives 
to implicit integration of even extremely stiff systems, and that these methods may permit 
integration of much larger networks than have been possible before in a number of fields. 



PACS numbers: 02.60.Lj, 02.30.Jr, 82.33.Vx, 47.40, 26.30.-k, 95.30.Lz, 47.70.-n, 82.20.-w, 
47.70.Pq 
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1. Introduction 

Problems from many fields of science and technology require the solution of large coupled 
reaction networks describing the flow of population between various sources and sinks. Some 
important examples include reaction networks in combustion chemistry Q, geochemical 
cycling of elements [4J, and thermonuclear reaction networks in astrophysics [ElO. The 
systems of differential equations that are commonly used to model these reaction networks 
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Figure 1. The CNO cycle. On the left side the main branch of the cycle is illustrated with 
solid arrows and a side branch is illustrated with dashed arrows. On the right side, the main 
branch of the CNO cycle is illustrated in more detail (/? stands for and a for ^He). 



usually exhibit stiffness, which we may think of loosely as arising from multiple timescales 
in the problem that differ by many orders of magnitude [i3l|71[8l|91. The most straightforward 
way to solve such a system might appear to be an explicit numerical integration scheme (for 
an explicit method, advancing a timestep requires only information known from previous 
timesteps). However, the textbooks routinely state [H [H |9l that such systems cannot be 
integrated efficiently using explicit forward finite-difference methods because of stability 
issues: for an explicit algorithm, the maximum stable timestep in a stiff system is typically set 
by the fastest timescales, even if those timescales are not of central interest. The traditional 
solution of the stiffness problem is to replace explicit integration with implicit integration 
(integration methods that require the values for derivatives at timesteps not yet evaluated). 

The astrophysical CNO cycle for conversion of hydrogen to helium (Fig. [I]) is an 
instructive example of the stiffness issue. In the CNO cycle the fastest rates typically are 
/3-decays with half-lives of order 100 seconds, but tracking main-sequence hydrogen burning 
may require integration of the hydrogen-burning network for as long as billions of years 
(~ 10^^ seconds). With explicit forward differencing the largest stable integration timestep is 
set by the fastest rates and will be of order 10^ seconds, so ~ 10^^ explicit integration steps 
could be required, even for the idealized case of constant temperature and density. Conversely, 
this same numerical integration requires at most a few hundred steps using implicit methods. 
For this reason, it is generally thought that explicit methods are not viable for extremely 
stiff networks. To quote Numerical Recipes [9J, "For stiff problems we must use an implicit 
method if we want to avoid having tiny stepsizes." 

Although implicit methods typically are stable for stiff systems, they require substantial 
additional computational overhead relative to explicit methods. For large coupled sets of 
equations, this normally takes the form of iterative solutions that require the inversion of large 
matrices at each step. The required matrix inversions dictate that, except where simplifications 
based on matrix structure can be exploited (for example, sparse-matrix methods), implicit 
algorithms may scale as poorly as quadratically or cubically with the size of the network. 
Thus, implicit methods can be expensive to implement for large networks, particularly those 
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that are coupled in real time to a broader problem such as hydrodynamical evolution. 

Despite the generally negative view of explicit methods for stiff systems sketched in 
the preceding paragraphs, they would be attractive options for complex networks if they 
could take larger timesteps because of their overall simplicity and highly-favorable linear 
scaling with network size relative to implicit methods. To increase the integration step size 
for explicit methods obviously requires overcoming formidable numerical stability problems 
that are associated with the stiffness. In principle, this might be accomplished by using 
approximate analytical solutions to reduce the stiffness of the equation set to be solved, 
thereby improving the overall stability of the network. In the first two papers of this series 
ID 13, we have shown in a variety of examples that this can be accomplished for systems 
that are not near equilibrium by using asymptotic and quasi-steady-state approximations to 
stabilize the numerical integration. However, as we shall now discuss, the nature of stiffness 
is different for networks that are near equilibrium. This implies that the modifications required 
to stabilize standard explicit methods when stiffness is encountered far from equilibrium must 
be altered when that same network approaches equilibrium. 

2. Varieties of Stiffness 

For the large and very stiff networks that are our primary interest here there are several 
fundamentally different sources of stiffness instability. Most textbook discussions of stiffness 
concentrate on an instability associated with small quantities that strictly should be non- 
negative being driven negative by an overly-ambitious numerical integration step. However, 
there are other guises that stiffness can take, as we now discuss. The equations that we must 
integrate take the general form 



where the yi{i =\ .. .N) describe the dependent variables (abundance in our examples), t is the 
independent variable (time in our examples), the fluxes between species i and j are denoted 
by {ff)i, and the sum for each variable i is over all variables j coupled to i by a non-zero 
flux {ff)i. For an N -species network there will be such equations in the populations j,, 
generally coupled to each other because of the dependence of the fluxes on the different yj. 

Because our discussion is intended to be quite general, we shall most often formulate the 
equations ^ in terms of generic population variables that are assumed to be proportional 
to the number density of species i. Where we give specific results for astrophysical networks, 
we shall use population variables most common for that field, such as the mass fraction X, and 
the (molar) abundance 7,, with 




= (/i+ - /r).- + (/2+ -f2),+---= £(/; - fi) 
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where A^a is Avogadro's number, p is the total mass density, A, is the atomic mass number 
and n, the number density for the species and by definition the mass fractions sum to unity 
if nucleon number is conserved: J^X, = 1. 

In Eq. ([U) the total flux has been decomposed into a component that increases the 
population and a component Fj that depletes the population yi, and in the third line this 
has been decomposed further into individual groups of terms — fJ. In the approach to 
equilibrium a species population entails a delicate balance between a total flux populating 
the species and a total flux depleting it. Near equilibrium the difference Fi = F-'^ — Ff^ can 
be orders of magnitude smaller than F-'^ or F- and small numerical errors in or F- can 
produce large errors in the difference Ft. Because of the population coupling in complex 
networks, this error propagates and compromises the accuracy of the network unless the 
timestep is short enough that the difference Fj is computed accurately for each population 
in the network. But this restriction means that the maximum timestep is set by the largest 
fluxes (that is, the largest stable timestep is determined by the inverses of the highest rates); 
this lands us back in the explicit integration conundrum that the maximum stable timestep is 
set by the fastest transitions, even if the primary interest is in quantities varying on a much 
longer timescale. 

Thus, in the approach to equilibrium the problem with explicit integration is not negative 
populations directly, but an unacceptable loss of accuracy that may occur even before any 
populations become negative. This is still a stiffness issue because it involves stability in 
systems with disparate timescales. In this case the contrasting timescales are the very rapid 
reactions driving the system to equilibrium relative to the very slow timescale associated 
with equilibrium itself. Thus any system nearing equilibrium can be expected to exhibit 
this form of stiffness instability. As we now address, this distinction among sources of 
stiffness is critical because these stiffness instabilities have fundamentally different causes 
and thus fundamentally different solutions. Furthermore, we shall find that the second class of 
instabilities can be divided into two subclasses requiring different stabilizing approximations. 
The approximations that we shall introduce in all these cases take care naturally of the 
first class of stiffness instabilities because they will prevent the occurrence of negative 
probabilities. 

3. Curing Stiffness in the Approach to Equilibrium 

We shall take equilibration to mean a situation where populations in the network are being 
strongly influenced by canceling terms on the right sides of the differential equations ([T]). In 
terms of the coupled set of differential equations describing the network, we may distinguish 
two qualitatively different conditions: 

(i) An equilibration acting at the level of individual differential equations that we shall call 
macroscopic equilibration. 

(ii) An equilibration acting at the level of subsets of terms within a given differential equation 
that we shall term microscopic equilibration. 
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Let us consider each of these cases in turn. 

3.1. Macroscopic Equilibration 

The differential equations that we must solve take the general form given in Eq. ([T]), dyi/dt = 
— Fj . One class of approximate solutions depends upon assuming that Fj^ ~ ~^ ^ 
(asymptotic approximations) or F^^ — Fj — > constant (quasi-steady-state approximations). 
We shall term this macroscopic equilibration, since these conditions involve the entire right 
side of a differential equation in Eq. ([U) tending to zero or a finite constant. In the first two 
papers of this series [1 , 2J we employed asymptotic and quasi-steady-state approximations that 
removed entire differential equations from the numerical integration for a network timestep by 
replacing them with algebraic approximate solutions for that timestep. These approximations 
integrate the full original set of differential equations, but they reduce the number of equations 
integrated numerically. This removal of equations from the numerical integration reduces the 
stiffness because it generally decreases the range of timescales in the numerical integration. 

3.2. Microscopic Equilibration 

In Eq. ©, Ff^ and F- for a given species / each consist of a number of terms depending on 
the various populations in the network, 

t='^."-'^r=E(/;-/7)- (3) 

At the more microscopic level, groups of individual terms on the right side of Eq. ([U) in 
the sum over j may come approximately into equilibrium (so that the sum of their fluxes 
is approximately zero), even if the macroscopic conditions for equilibration are not satisfied 
and asymptotic or quasi-steady-state approximations are not well justified for the species i. 
The simplest possibility is that forward-reverse reaction pairs such as A -|- 5 C, which 
will contribute flux terms with opposite signs on the right sides of differential equations 
in which they participate, come approximately into equilibrium. As we shall demonstrate, 
this introduces new (often fast) timescales into the problem that are at best only partially 
removed by asymptotic and quasi-steady-state (QSS) approximations. The new sources of 
stiffness associated with these microscopic equilibration effects explain why asymptotic and 
QSS approximations remain accurate but their timestepping becomes very inefficient in the 
approach to equilibrium. 

As we elaborate in the remainder of this paper, this new source of stiffness associated 
with close approach to microscopic equilibration requires a new algebraic approximation 
that removes groups of such terms from the numerical integration by replacing their sum 
of fluxes with zero. Such an approximation will not generally reduce the number of equations 
to be integrated numerically in a network timestep, but can (dramatically, as we shall see) 
reduce the stiffness of those equations by systematically removing terms with fast rates from 
the equations. This reduces the disparity between the fastest and slowest timescales in the 
system. Thus shall we convert the approach to equilibrium in an explicit integration from a 
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liability into an asset. As part of this elaboration we shall also reach two important general 
conclusions: (1) Approximations based on microscopic equilibration are much more efficient 
at removing stiffness than those based on macroscopic equilibration, because they more 
precisely target the sources of stiffness in the network. (2) Macroscopic and microscopic 
approximations can complement each other in removing stiffness from the equations to be 
integrated numerically and thus the two together are more powerful than either used alone. 

4. Methods for Partial Equilibrium 

Let us begin to develop some methods to deal effectively with explicit integration in the 
approach to equilibrium. In doing so we draw substantially on the work of David Mott [fTOll . 
but we will extend these methods and obtain much more favorable results for extremely stiff 
networks than those obtained in the original work of Mott and collaborators. 

The basic idea of partial equilibrium (PE) methods is to inspect the source terms f^^ 
and /-^ associated with individual reaction pairs in the network for approach to equilibrium 
(instead of the composite flux terms F-'^ and that are the basis for asymptotic and quasi- 
steady-state approximations). Once a fast reaction pair nears equilibrium, its source terms 
are removed from the direct numerical integration in the ordinary differential equations and 
its effect is incorporated through an algebraic constraint implied by the equilibrium condition. 
Those reactions not in equilibrium still contribute to the net fluxes for the numerical integrator, 
but once the fast reactions associated with an equilibrating reaction pair are decoupled 
from the numerical integration the remaining system typically becomes much less stiff. To 
illustrate, consider a representative 2-body reaction, 

a + b^c + d. (4) 

The source term for this reaction pair can be expressed in the general form 

fabled = ^{fa+h^c+d — fc+d^a+b) = ±(%a>'fc " ^rjcjd), (5) 

where the yi denote population variables for the species i and the ks are rate parameters. This 
source term will contribute to the right side of the differential equations describing the change 
in population for all four species a, b, c, d: 

fabled + Other terms changing ya (6) 
fabled + other terms changing yb (7) 

- fabled + Other terms changing y^ (8) 

- fabled + other terms changing yj (9) 

Thus removing or reducing the stiffness associated with this single reaction pair influences a 
whole set of populations and associated reactions, and so can reduce the overall stiffness of 
the system. Strictly, the reaction pair of Eq. (U) is in equilibrium if fabled = 0- Of greater 
interest will be partially-equilibrated systems, where some reactions maintain / ~ as the 



dya 
dt 
dyb 
dt 
dye 
dt 
dyd 
dt 
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system evolves but others have f j^O (and the members of these two sets may change over 
time). Our goal will be to develop methods to integrate such systems using asymptotic or 
quasi-steady-state algorithms, but for modified equations of reduced stiffness in which the 
flux sums for the / ~ reactions have been replaced by equilibrium constraints. 

From Eqs. ([l])-® the single reaction pair a + b c + d appears to have four 
characteristic timescales associated with the rate of change for the four populations ya, yb, 
yc, and y^, respectively. For a + b ^ c + d considered in isolation, the differential equation 
governing the abundance of species a is 

^ = F+ - F- = Ky^y, - hy^y^ = - ^ , 

dt Ta 

where Ta = ^/(kfyb)- If we assume the abundances of c and d and the rate parameters to 
remain approximately constant in a timestep, the timescale Ta characterizes the rate of change 
of the species a. Likewise, for the other three components of the reaction pair we may write 
similar differential equations and define similar timescales, 

1111 

% = T^c = T^d 



hjb hja kyd Kyc 

Our first task is to construct a single timescale that characterizes equilibration of a reaction 
pair such as a+b ^ c + d, considered in isolation. To do so we introduce the idea of conserved 
scalars [llOjl. 



4.1. Conserved Scalars and Progress Variables 

As a simple initial illustration, consider the reaction pair a # 2b, which has a source term 

fa^2b = hya-kryl- (10) 

If no other reactions altered the populations ya and yb, 
dt dt 

Thus Idyajdt -Vdybldt = and lya+yb = constant. The quantity lya-Vyb is an example of 
a conserved scalar. It is conserved by virtue of the structure of a lb, not by any particular 
dynamical assumptions; thus conservation of this quantity is independent of whether the 
reaction is near equilibrium or not. It is convenient to introduce new variables that are the 
difference between the initial values of yi and their current values 

^ya = ya-y^a 8yb = yb-yl- 

The initial values j° and y^ are constants so the differential equations (fTT)) can then be written 



as 



dSya dSyb 
dt dt 
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This suggests defining a progress variable A for the reaction characterized by fa^2b which 
satisfies 

-^=fa^2b Ao = A(r = 0)=0. (13) 
By comparing Eq. (fT3l) with Eq. (fT2l) we have A = — = j5yt„ so that 

= Sya = -A +3;° jfc = Sjfc = 2A +3;^. (14) 

Thus Eqs. (flOl) and (fTTI) may be rewritten in terms of the single variable X and the solution of 
two differential equations for two unknowns is reduced to the solution of Eq. (fT3l) for a single 
unknown A [which may then be used to compute ya and yi, through Eq. (fT4l)1. 

Let's now apply these ideas to the general 2-body reaction a + b^c + d of Eq. (U). 
Assuming conservation of particle number, it is clear that the following constraints apply to 
this reaction 

ya-yb = ci ya+yc = c2 ya+yd = c3, (15) 

where the are constants. Losing one a by this reaction requires the simultaneous loss of 
one b, so their difference must be constant, and every loss of one a produces one c and one 
d, which explains the second and third equations. These constraints follow entirely from the 
structure of the reaction and are independent of dynamics. The constants can be evaluated by 
substituting the initial abundances into Eq. (flSl) . 

c,=y'a-yl C2=y'^+y^. c,=yl^yl 
The differential equation for ya is 

^ = -hyayb + hycyd , (16) 

dt 

which can be rewritten using Eq. (fT5l) as 

^ = ayl+bya + c, (17) 

where 

a = k^ — kf b = —k^{c2 + C3) + fcfCi c = k^C2C3 . 

As we demonstrate below, the approach to equilibrium for any 2-body reaction pair can be 
described by a differential equation of this form, and the approach to equilibrium for any 3- 
body reaction pair can be approximated by a differential equation of this form. If a = the 
solution of this differential equation gives the quasi- steady- state (QSS) solution, which we 
have already examined in the second paper of this series [2]. For the general case a ^ 0, let 
us define a quantity 

q = 4ac-b^ <0. (18) 

The case ^ = is the trivial solution where all abundances are zero, so we are interested in 
solving Eq. (fTTI) for negative values of q. The general solution for a 7^ and ^ < is [fTOl[TT| 

1 /, , l + 0exp(-v/^O \ 

ya{t) = -:^(b + V^ I ^ , (19) 

2a V I — ^exp{—y/^t) J 




Figure 2. (a) Time evolution of the solution ( fT9] l assuming a, b, and c to remain constant. The 
characteristic timescale for approach to equilibrium defined by Eq. ( |23] | is labeled T and the 
equilibrium value of y{t) defined by Eq. ( 1211 1 is denoted by y. To illustrate we have assumed 
the initial value = 0. For times considerably larger than T the general solution ( fT9] l saturates 
at the equilibrium solution ( 1211 1. (b) Behavior of the exponential factor in Eq. ( fT9l l. 



where 

The equilibrium solution corresponds to the limit ? — > oo of Eq. (fT9l ): 

3J,^3;^q = -l(Z7 + v^). (21) 

By analogy with the earlier discussion we define a progress variable 

l{t) = -ya{t)+yl. (22) 

Once ya has been determined the constraints (fT5l) may be used to determine the other 
abundances: 

yb{.t)=ya{t)-ci yc{t) = c2-ya{t) yd{t) = c3-ya{t)- 

The approach of the reaction a + b ^ c + d to equilibrium is thus governed by a single 
differential equation ([IT]) , which may be expressed in terms of either a single one of the 
abundances yi, or the progress variable A defined in Eq. ((22)) . The general solution of 
this differential equation is given by Eq. (fT9l ). from which the rate at which the reaction 
a + b ^ c + d evolves toward the equilibrium solution (121]) is determined by a single timescale 

^ (23) 



These equilibrium timescales are illustrated in Fig. [2l We may then estimate whether a 
reaction is near equilibrium at time t by requiring 



yi 



< £: (24) 
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for each species / involved in the reaction, where is the actual abundance, j,- is the 
equilibrium abundance computed from Eq. (|2T]) . and the user-specified tolerance can 
depend on i but will be taken to be the same for all species in the simplest implementation. 
Alternatively, the equilibrium timescale (|23l) in comparison with the current numerical 
timestep may be used as a measure of whether a reaction is near equilibrium. If T is much less 
than the timestep, it is likely that equilibrium can be established and maintained in successive 
timesteps, even if it is being continually disturbed by other non-equilibrated processes. 

4.2. Reaction Vectors 

In a large reaction network there could be thousands of reactions to be examined for their 
equilibrium status at each timestep, so implementation of a partial equilibrium approximation 
requires a substantial amount of bookkeeping. Our task will be facilitated by systematic ways 
to catalog and examine the equilibrium status of reactions. We employ a formalism, adapted 
from the thesis work of David Mott [fTOll . that exploits the analogy of a reaction network to a 
linear vector space. This will have two large advantages for us: (1) It will place at our disposal 
well-established mathematical tools. (2) Treating the reaction network as a linear vector space 
will permit formulation of a partial equilibrium algorithm that is not tied closely to the details 
of a particular problem, thus aiding portability across disciplines. 

It is useful to view the concentration variables for the n species A, in a network as 
components of a composition vector 

y = {yi,y2,y3,---yn), (25) 

which lies in an n-dimensional vector space 4>. The components can be any parameters 
proportional to number densities for the species labeled by A,, and a specific vector in this 
space defines a particular composition. Any reaction in the network can then be written in the 
form 

£aA-^£M/ (26) 

i=l i=l 

for some sets of coefficients {a,} and {bt}. For example, consider the CNO cycle of Fig. [T] 
and choose an ordering 

y = (p = iR, a = ^He, 12c, 13n, ^^C, ^^N, ^^O, ^^N) . (27) 

(Note that usually the positrons /3+ and the gamma rays 7 emitted in Fig. [T] are not tracked 
explicitly in the network.) Letting the composition variables correspond to the mass fraction 
Xi for a species, the ^^C{p, yY^N reaction (that is, ^11+ ^^C ^^N) corresponds to Eq. (|26|) 
with a = {l, 0,1, 0,0,0,0,0} and Z7 = {0,0,0,1, 0,0,0,0}. The coefficients on the two sides of 
a reaction may be used to define a vector r e 4> that has the form 

r = {bi-aiM-ai^-'-^bn-an), 

and defines a composition displacement vector in 4> associated with the reaction. For example, 
for the reaction -f _i. 13]^ jj^g components of r are (- 1 , 0, - 1 , 1 , 0, 0, 0, 0) . 



Explicit Integration of Extremely-Stiff Reaction Networks: Partial Equilibrium Methods 1 1 



4.3. Conservation Laws 

Given an initial composition yo = (jj, J2'-y3' • ■ composition that can be produced by 

a single pair of reactions labeled by i is of the form y = yo + oCiVj, where a, is some scalar 
quantity and for a set of k reactions the final composition must be of the form 

k 

y = yo + Y, ^i'^' (28) 

Define a time-independent vector c = (ci , C2, C3, . . . c„) G 4> that is orthogonal to each of the k 
vectors in Eq. (|28]) . From Eq. ([281) 

k 

c- (y-yo) = c- J^a^r, = 0, 

which implies that c - y = c - yo and therefore that 

n n 

CiJi = '^iyo = constant. (29) 

Thus any vector c orthogonal to the reaction vectors r\, r2, ■ ■ .r^ gives a linear combination 
of species abundances that is invariant under the reactions defined by the vectors ri. Equation 
(|29l ) defines a conservation law following only from the structure of the network, independent 
of any dynamical considerations, and thus must be valid irrespective of dynamical conditions 
in the network. The conserved quantities may be determined by forming a matrix r having 
rows corresponding to the reaction vectors r,, and solving the matrix equation r ■ c = for the 
vector c. 

4.4. Example: CNO Cycle 

Let us illustrate some of the preceding ideas by applying the classification scheme to the 
network describing evolution of the main part of the CNO cycle with time. From Fig.[Tl the 
reactions of the closed cycle are 

1: p + ^^C i^N + y 2: ^ i3c + /3+ + Ve 

3: p + i^C ^ i^N + y 4: p^-^^^ ^ + y (30) 

5: 1^0 ^ + + 6: p + ^^N ^ a + ^^C 

which we shall assume to not be reversible under the temperature and density conditions 
characteristic of the CNO cycle. (Since the reactions are assumed to not be reversible, we 
shall not apply the partial equilibrium approximation to this set of reactions. But that is 
irrelevant for the classification scheme, which is independent of dynamical conditions in the 
network.) Utilizing the vector (|27l) . the source terms for the reactions (l30l) are 

1 : Fi3 = ^i3j 2 : F4 = ^4/ 3 : F15 = ^15^^^ 
4 : F16 = ki(,y^y^ 5 : F7 = ^7/ 6 : Fig = ki^y^y^. 
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Table 1. Reaction classes in the REACLIB |[T2| library 



Class Reaction Description or example 



1 


a ^ b 




/3 -decay ore capture 


2 


a ^ b + c 


Photodisintegration + a 


3 


a — > b + c + d 


12c ^ 3a 


4 


a + b — 


> c 


Capture reactions 


5 


a + b — 


> c + d 


Exchange reactions 


6 


a + b — 


> c + d + e 


2h + ^Be ^ iH + 2^He 


7 


a + b — 


>-c + d + e + f 


3He + ^Be^2iH + 2'*He 


8 


a + b + c-^d (+e) 


Effective 3-body reactions 



where the indices refer to the positions in the vector (|27]) and the ^s are rate parameters that 
in the general case depend on temperature and density. With the ordering (1271 ) the reaction 
vectors corresponding to Eq. (|30l) are 

n = (-1,0, -1, 1,0,0,0,0) r2 = (0,0,0, -1, 1,0,0,0) 

r3 = (-1,0,0,0,-1, 1,0,0) r4 = (-1,0,0,0,0,-1, 1,0) (31) 

r5 = (0,0,0,0,0,0,-1,1) r6 = (-1,1, 1,0,0,0,0,-1) 

By forming a matrix r having rows given by the reaction vectors r,- and solving the matrix 
equation r ■ c = for the vector c by gaussian elimination, we find two useful conservation 
laws, specified by the vectors ci and C2. The first corresponds to conserving the total number 
of neutrons plus protons (conservation of nucleon number). The second corresponds to 
conservation of the sum of number densities for all the carbon, nitrogen, and oxygen (CNO) 
isotopes in the network, which represents an elegant derivation of the well-known property 
that the CNO isotopes catalyze the CNO cycle, but are not consumed by it. 

4.5. Reaction Group Classes 

It will prove useful to associate inverse reaction pairs in reaction group classes (reaction 
groups, or RG for short). We employ the individual reaction classifications used in the 
REACLIB library fl2l| that are illustrated in Table [TJ In this classification reactions are 
assigned to eight categories, depending on the number of nuclear species on the left and right 
side of the reaction equation. From this classification it is clear that there are five independent 
ways that the reactions of Table [U can be combined to give reversible reaction pairs. This 
forms the basis for the reaction group classification illustrated in Table [2l For example, 
reaction group class B consists of reactions from REACLIB reaction class 2 (a b -i- c) 
paired with inverse reactions (b -i- c — t- a), which corresponds to REACLIB reaction class 4. 

Reaction group class A corresponds mostly to /3-decays and their inverses, so it is 
important in partial equilibrium calculations only if neutrino reactions are included. There 
are only a few reaction pairs of broad importance in classes C and E other than triple- 
a (3a ^^C), so in many practical applications the most important reaction groups lie 
predominantly in reaction group classes B and D. The reason that REACLIB class 8 appears 
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Table 2. Reaction group classes 



Class 


Reaction pair 


REACLIB class pairing 


A 


a^b 


1 with 1 


B 


a + b ^ c 


2 with 4 


C 


a +b + c ^ d 


3 with part of 8 


D 


a + b ^ c + d 


5 with 5 


E 


a + b^c + d + e 


6 with part of 8 



in both reaction group classes C and E is the ambiguity in its definition in Table [H since it 
contains reactions that can have either one or two nuclear products on the right side. Notice 
that reaction group class D is composed of REACLIB class 5 paired with itself, and that 
there are no reaction group classes that involve REACLIB reaction class 7 because it has four 
nuclear species on the right side and REACLIB contains no reactions with four bodies on the 
left side of the reaction equation to pair with it. 

Rare exceptions to these observations can occur if a reaction involves a catalyst (same 
species appearing on both sides of the equation with the same coefficients). For example, 
consider the reactions n + 2^He # ^Be and p + ^Be n + p + 2^He. The first reaction pair is 
classified as reaction group class C since it pairs a REACLIB class 8 reaction with a REACLIB 
class 3 reaction. The components of this reaction pair have reaction vectors that differ from 
each other by a sign. The p + ^Be reaction is REACLIB class 7, which normally does not 
have a REACLIB class to pair with in equilibrium. However, the net effect of this reaction on 
populations is the same as ^Be — )■ n + 2'*He since the proton appearing on both sides of the 
equation is catalytic and cancels in the population changes. Thus p + ^Be n + p + 2'*He has 
the same reaction vector as ^Be — t- n + 2^He and is effectively also the inverse of n + 2^He. The 
only other such anomalous example contained in REACLIB [12] is afforded by the reactions 
n + p ^ d and p + d ^ n + p + p, where d denotes the deuteron ^H. 

For each reaction group class the differential equation governing the reaction pair takes 
the form given by Eq. (flTl) . dy/dt — ay^ + by + c, where y is either a variable proportional to a 
number density for one of the reaction species, or a progress variable that measures the change 
in initial abundances associated with the reaction pair, and the coefficients a, b, and c will be 
assumed constant within a single network timestep. An exception occurs for reaction group 
classes C and E, which contain 3-body reactions so that the general form of the differential 
equation involves cubic terms dy/dt = ay^ + /3y^ + 7> + e . One could take a similar approach 
as before, solving this cubic equation for the partial equilibrium properties for the reaction 
group. However, these "3-body" reactions in astrophysics are typically actually sequential 
2-body reactions and we employ an approximation that in any timestep y{t)^ ~ y'^^^y{tY, 
where the constant the value of y{t) at the beginning of the timestep. This reduces 

the cubic equation to an effective quadratic equation of the form (fTTl) . with a = ay^*^) + /3, 
b = Y, and c = e. For the case a 7^ we have already described the corresponding general 
solution, associated timescale, equilibrium abundances, and test for approach to equilibrium 
of the reaction in Eqs. (fr7l) - (|24)) . Our tests suggest that this is a very good approximation in 
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typical astrophysical environments and we shall treat all 3-body reactions as effective 2-body 
reactions. 



4.6. Equilibrium Constraints 

If a reaction pair from a specific reaction group class of the form (|26|) is near equilibrium, 
there will be a corresponding equilibrium constraint [fTSlI 

n?''""'^=if, (32) 

where K is some ratio of rate parameters. For example, consider the reaction group class E 
pair a + b"^ c-\-d-\-e, with differential equations for the populations ji 

At equilibrium, the requirement that the forward flux —kfyayb and backward flux ^rjcjdje in 
the reaction pair sum to zero implies the constraint 

yc-ydye h 

which is of the form (|32l) . 



4. 7. Reaction Group Classification 

Applying the principles discussed in the preceding paragraphs to the reaction group classes 
in Table |2] gives the results summarized for reaction group classes A-E in [Appendix A 



This gives us a complete classification scheme as a basis for applying a partial equilibrium 
approximation to arbitrary astrophysical thermonuclear networks. However, the methodology 
is of broader significance. First, since any reaction compilation in astrophysics could 
be reparameterized in the REACLIB format, this classification scheme provides a partial 
equilibrium bookkeeping for any problem in astrophysics. Second, for any large reaction 
network in any field, the classification techniques illustrated here can be applied to group all 
reactions into reaction group classes, and to deduce for each reaction group class the quantities 
necessary for applying a partial equilibrium approximation. All that is required is to formulate 
the network as a linear algebra problem by choosing a set of basis vectors corresponding to 
the species of the network, and then to define the corresponding reaction vectors within this 
space. (Mathematically the choice of the vector space is arbitrary as long as it provides a 
faithful mapping of possible species and reactions, but physical interpretation may be aided 
by judicious choices in specific fields.) In principle this need only be done once for the 
networks of importance in any particular discipline. 



5. Example Illustrating the Basic Idea 

Let us work through a simple example illustrating why partial equilibrium, and partial 
equilibrium in conjunction with the explicit asymptotic method, could greatly reduce the 
stiffness associated with numerical integration of a set of coupled differential equations. 
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Table 3. Reaction groups and constraints 



Group Reactions Class Constraint 



1 


fl4 




c 


B 


yayt 


4" 
" 4" 


2 


3a 


^b 




C 


yl = 


jriyh 


3 


a 4 


c ^ 


d 


B 


yayd 





5.1. Network and Reactions 

For this example we consider the simple network a ^ b ^ c ^ d, including the reaction 
groups 

a + b^c 3a ^b a + c^d. (33) 

Our general conclusions will not be specific to astrophysics, but in fact these equations have 
the form of the first part of an astrophysical thermonuclear alpha-particle network 

with oc-capture ( a + b ^ c), photodisintegration (c — j- a + b), and the triple-a reaction 
(3a ^ ^^C) included. From the reaction group classification in [Appendix A[ these reaction 



groups belong to classes B and C, and the associated equilibrium constraints are given in 
Table [3j where the k are the rate parameters (generally time-dependent), with the superscripts 
denoting the reaction group and the subscripts indicating forward (f) and reverse (r) directions 
in Eq. (1331) . Thus ki^^ is the rate for c ^ a + b. The coupled set of ordinary differential 
equations to be solved is 

_ ,(2) 3 ,,(2)^, k^^K, ^, ^k^^K, ,(3) 

— —kf + % — k^ yaJb + Jc — Jajc + /^r K^V 

— '^f ya~ '^i yh + % yc - yayt ( j j) 

,(i) J3)„,, nf.\ 
-j^ — yayt — % yc — yayc + % yd (jo) 

dyd _ , (3)^, , (3) 

-7- — yayc-ki yd- (-?/) 

dt 

From Table [3] there are potentially three constraints available if the system comes into full 
equilibrium, and one or two constraints if it is in partial equilibrium. In addition to the 
equilibrium constraints, we may wish to impose constraints associated with conservation laws 
for the system, such as preservation of particle number. In actual applications this will be very 
important, but in this example we ignore conservation laws and concentrate on understanding 
the role of equilibrium constraints. There are two general approaches that we might take 
to using equilibrium constraints to simplify the solution of the differential equations in Eqs. 
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(i) Use the constraints to reduce the number of equations that we have to solve in a given 
timestep. 

(ii) Use the constraints to reduce the stiffness of the equations that we have to solve in a 
given timestep. 

The first approach reduces the number of equations to solve, but in most cases does not change 
the stiffness much for the equations that are solved until substantial numbers of equations 
have been removed from the numerical integration. In the second approach, we still solve the 
same number of equations, but the equations that we solve are less stiff than the original ones 
because we have modified the structure on the right side of the equalities in Eqs. (I34l)-(l37l). 
We shall see that the second approach tends to naturally remove the fastest timescales that 
remain in the system when it is applied, so it can have a dramatic effect on the stiffness of 
the system. Thus in this paper we shall only outline using constraints to reduce the number of 
equations, and then concentrate on how constraints can be used to reduce stiffness. 



5.2. Using Equilibrium Constraints to Reduce the Number of Equations 

Consider the constraint from a + b ¥^ c written in the form jayb — Kjc = 0, where 
k^"^^ /k^^\ Taking the derivative with respect to time of this expression gives 

djb ^ 

dt dt dt dt ' 



which is a constraint that could be used to eliminate one of Eqs. (|34l) - (|37l) from the numerical 
integration, say Eq. (|35l) . Notice that removing Eq. (l35l) leaves all of the rate parameters of the 
original problem present in the remaining equations, so it presumably has had only a small 
impact on the stiffness. In a similar manner, as other reaction pairs come into equilibrium 
the associated constraints can be used to remove additional equations from the numerical 
integration. We shall not use this approach in the present context, but we note in passing 
that this represents a systematic way to introduce a partial equilibrium approximation into an 
implicit-method calculation. In that case, reducing the number of equations to be integrated 
can be significant because it reduces the sizes of the matrices that must be inverted at each 
integration step. 



5.3. Using Equilibrium Constraints to Reduce Stiffness 

Instead of using the equilibrium constraints to reduce the number of equations, let us now 
outline how to use them to reduce the stiffness of the original set of equations by applying the 
constraint directly to the abundances rather than their derivatives. Suppose that the reaction 
a-\-b^ c comes into equilibrium, implying the constraint jaj^ = k[^^yc/k^^\ Substituting this 
for yajb in Eqs. (I34l)-(l37l). various terms cancel and we obtain the modified equations 



^=-4^V2+tr-v»-*fw+tf'y. (38) 

'-^^kfhl-k^b. (39) 
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!|=-*fV* + tf'>. (40) 

^|l = tfW-t<^,. (41) 

This is the same number of equations as in (I34l)-(l37l). but now k^^^ and k[^^ no longer appear 
on the right sides of the differential equations. Furthermore, we may expect that the rate 
parameters that have been removed were fast ones, because they were responsible for bringing 
the first reaction group into equilibrium in the network. Therefore, on general grounds we may 
expect that the differential equations Eqs. (l38l)- (|4T]) are less stiff than the original equations 
(I34l)-(l37l). because there will be less disparity between the fastest and slowest rates in the 
network. 

Let us suppose further that when the fluxes on the right side of Eqs. (I38l)-(|4T1) are 
computed we find that two of the differential equations — let's say (|38l) and (|39l ) — satisfy the 
asymptotic condition. Then these two entire equations would be solved for the timestep by 
the algebraic asymptotic formulas, which eliminates another whole set of terms involving 
different rate parameters from the the numerical integration. Thus we can see conceptually 
how simultaneous use of partial equilibrium and asymptotic approximations could reduce 
considerably the effective stiffness of a complex set of equations. Continuing our example, 
suppose that the reactions a-\-b^ c remain in equilibrium (if a reaction group in equilibrium 
drops out of equilibrium at a later time, the corresponding flux terms must be restored to the 
differential equations) and at a later time the reaction 3a ^ b also comes into equilibrium. 
Now there are two constraints to be applied to Eqs. (I34l)-(l37l). or one new one to be applied to 
Eqs. (l38l)- (|4T]) . Applying the new equilibrium constraint yl = k^^^h/kf^ to Eqs. (|38l) - (|4T]) . 
we obtain 

^-^=-kf\ayc + k?\, {AD 

^ = (43) 
dt 

+ (44) 

5|^ = *fW-*PV.. (45) 
Again we have the same number of equations to integrate numerically as originally (although 
Eq. (|43]) has become trivial), but additional fast terms have been removed from the right side 
of the differential equations, reducing their stiffness even further. As before, if any equations 
in the preceding set satisfy the asymptotic condition, these equations may be removed from 
the numerical integration in favor of an algebraic solution, potentially further reducing the 
stiffness of the numerical system. 

Finally, suppose that a + b ^ c and 3a ^ b remain in equilibrium, and at some later time 

(3) (3) 

a + c^ d comes into equilibrium too, implying a new constraint yayc = k) yd /kf . Inserting 



this expression for yayb into Eqs. (|42l) -(l45l). we obtain the system characteristic of complete 
equilibrium, 

dya dyh dye dyd ^ 
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Of course this set of equations can be integrated trivially, but that is the point! The 
systematic application of partial equilibrium techniques to an intrinsically highly- stiff system 
has produced an approximately equivalent system in which all numerical stiffness has been 
removed (since no timescales remain in the equations). Thus there is no stability restriction 
on the timestep, should we choose to integrate Eq. (l46l) numerically. 

In realistic large and stiff networks we will generally be integrating numerically in 
regimes where at least some reactions are not fully in equilibrium, so the trivial limit 
of Eq. (|46l) will seldom be reached except for systems very near overall equilibrium. 
Nevertheless, the preceding example illustrates that large reductions in stiffness may still be 
realized through systematic application of the constraints for those reactions that do come 
into equilibrium. In realistic networks we may expect that in the approach to equilibrium, 
situations where the numerical system being integrating is only somewhat removed from a 
trivial one of the form given by Eq. (l46l) can occur frequently. In those cases, we may expect 
large increases in the maximum explicit integration timestep from systematic application of 
the constraints implied by reactions that come into equilibrium, coupled with asymptotic 
or QSS approximations employed when conditions warrant it, to the resulting simplified 
equations. 

6. General Methods for Partial Equilibrium Calculations 

We now have a set of tools to implement partial equilibrium approximations, but there are a 
number of practical issues that require resolution before we can make realistic calculations. 
To that end, let us now outline a specific approach to applying partial equilibrium methods. 

6.1. Overview of Approach 

The partial equilibrium method will be used in conjunction with the asymptotic approximation 
in the following examples. (Although we shall not address it in this paper, a similar algorithm 
could be employed that replaced the asymptotic approximation with the quasi- steady- state 
approximation described in Ref. [2J.) That is, we shall remove fluxes from the numerical 
integration corresponding to reaction pairs that are in equilibrium, but the integration is still 
performed formally over the full reaction network, subject to an asymptotic approximation. 
Once the reactions of the network are classified into reaction groups, the algorithm has three 
basic steps: 

(i) During a numerical integration step one begins with the full network of differential 
equations, but in computing the net fluxes all terms involving reaction groups that are 
judged to be equilibrated (based on criteria determined by isotopic populations at the 
end of the previous timestep) are assumed to sum identically to zero net flux and are 
omitted from the flux summations. 

(ii) A timestep At is then chosen (using a variant of the timestepping algorithm described 
in Refs. [fT] l2l|), and this is used in conjunction with the fluxes to determine how many 
isotopes in the network satisfy the asymptotic condition according to the criteria of Ref. 
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[[T]|. For those isotopes that are not asymptotic, the change in abundance for the timestep 
is then computed by ordinary forward (explicit) finite difference, but for those isotopes 
judged to be asymptotic the abundance changes for the timestep are instead computed 
using the analytical asymptotic approximation of Ref. [Jj. 

(iii) Finally, for all isotopes in reaction groups taken to be in equilibrium at the beginning of 
the current timestep, it is assumed that reactions not in equilibrium will have driven these 
populations slightly away from their equilibrium values during the numerical timestep. 
These populations are then adjusted, subject to a condition that the sum of the mass 
fractions remain equal to one, to restore their equilibrium values at the end of the 
timestep. 

Hence the partial equilibrium approximation does not reduce the number of equations to be 
integrated but instead removes the stiffest parts of their fluxes in each timestep. In contrast, the 
asymptotic approximation reduces the number of differential equations integrated numerically 
within a timestep by replacing the numerical forward difference with an analytically- 
computed abundance for those isotopes satisfying the asymptotic condition. 

6.1.1. Reducing Stiffness Both of these approximations reduce stiffness in the integration by 
replacing numerical finite difference with algebraic conditions, but by different and potentially 
complementary means. The partial equilibrium method operates microscopically to remove 
individual components of reaction fluxes that become stiff in the approach to equilibrium; the 
asymptotic approximation operates macroscopically to remove the entire net flux altering the 
abundance of individual isotopes from the numerical integration. The partial equilibrium and 
asymptotic approaches are potentially complementary because partial equilibrium can make 
the right side of the differential equation for a given isotope less stiff, even if the isotope does 
not satisfy the asymptotic condition, while conversely the asymptotic condition (by removing 
the entire flux of selected differential equations from the numerical update) can effectively 
remove stiff reaction components even if they do not satisfy partial equilibrium conditions. 
For compactness we shall often refer simply to the partial equilibrium (PE) approximation in 
what follows, but it should be understood that this means partial equilibrium approximation 
plus asymptotic approximation as described above. 

6.1.2. Operator-Split Restoration of Equilibrium The third step (restoration of equilibrium) 
in the above algorithm may be viewed as an operator-split separation of physical timescales 
within a single numerical network timestep (with evolution of the network itself already being 
treated on a different level as operator- split from the evolution of the hydrodynamics). The 
idea is illustrated schematically in Fig.[3l Establishment of equilibrium for individual reaction 
groups assumed to be in equilibrium within a single numerical timestep is fast compared 
with the timescale for reactions causing net changes in the abundances within a timestep. 
Thus, if we were to do the integration exactly, equilibrium for those reaction groups would 
be maintained during the timestep; but then we would have a highly-stiff numerical system 
combining fast and slow components. Our approximation replaces this actual evolution by a 
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Yn = Yn yn+^ 




Operator-split timestep Af 



Figure 3. Illustration of the "operator-split" restoration of equilibrium abundance for an 
isotope. At the beginning of a network timestep that advances from /„ to the abundance yn 
is assumed to take its equilibrium value y,, . An explicit numerical timestep Af = f„+i — f„ is then 
taken for all abundances, but with the contribution of equilibrated reaction groups removed 
from the fluxes. The isotope in question starts the timestep with its equilibrium abundance, 
but the numerical timestep will generally leave the abundance at a value y„+\ that is not equal 
to the new equilibrium value y„^\ calculated from Eq. UaS because of its coupling to non- 
equilibrium reactions (the effect of the diagonal arrow). If the equilibrium approximation is to 
be maintained for the next timestep, the equilibrium value y-n+i must be restored before that 
timestep is taken. We could in principle restore the equilibrium by integrating numerically in 
a second step (the dashed blue line). However, we already know the result of that integration: 
it is given by Eq. (I2TI 1. because by hypothesis Af is larger than the characteristic equilibration 
time for all reactions assumed to be in equilibrium. Thus we accomplish the second operator- 
split step algebraically rather than numerically by replacing the numerically-computed y„+i 
with y„+i specified through Eq. (ISTT l (the vertical arrow on the right side). 

two-Step process: 

(i) Evolve all network abundances by explicit forward differencing and asymptotic 
approximations, with the equilibrium reactions assumed to maintain their equilibrium 
(net zero) fluxes over the timestep. (This is what justifies removing the fluxes that are 
assumed equilibrated from the flux summation.) 

(ii) Then, at the end of this first step, evolve all abundances assumed to be participating in 
equilibrium (which will have been disturbed by non-equilibrium processes in the first 
step) from their values computed at the end of the first step back to equilibrium values 
using Eq. (|2TI) . while holding non-equilibrium abundances constant. 

Unlike the case for the usual operator- splitting between evolution of the hydrodynamics and 
evolution of the reaction network, this second step is not computed numerically by finite 
difference but is instead implemented through algebraic constraints. This is possible because 
the equilibrium assumption for a given reaction group means that it should evolve to the 
equilibrium solution before the end of the timestep. Therefore, since we know how the story 
ends, rather than integrating it explicitly we may simply replace the population at the end of 
the numerical timestep by its equilibrium value, calculated from Eq. (|2T)) . We gain from this 



Explicit Integration of Extremely-Stiff Reaction Networks: Partial Equilibrium Methods 2 1 

separation of timescales and solution of the fast timescale by algebraic constraint rather than 
a finite-difference approximation a potentially large reduction of stiffness for the remaining 
part of the system that is integrated numerically. We shall demonstrate below that this can 
increase by orders of magnitude the maximum stable and accurate timestep for the overall 
integration in near-equilibrium conditions. 

6.1.3. Complications in Realistic Networks The complication for this basic idea in a realistic 
network with more than one reaction group in equilibrium is that there may be more than one 
computed equilibrium value for a given species that is assumed to be in equilibrium. This is 
because the equilibrium abundance Eq. (|2T1) must be computed separately for each reaction 
group, and an isotope will generally be found in more than one reaction group. For example, 
assume an a network in which both reaction group I, defined by a + ^^Si ^^S, and reaction 
group II, defined by a+^^Ca. ^ ^Ti, are assumed to be in equilibrium. Then by our algorithm 
Ya should assume its equilibrium value at the end of each timestep, but for each timestep there 
will be two values for the equilibrium abundance Ya, corresponding to Eq. (|2T1) solved for 
reaction group I and reaction group II, respectively. In the general case these values could 
be different, since the rates entering into Eq. (|2TI) are different for the two reaction groups, 
and equilibrium within each reaction group is specified only to the tolerance implied by in 
Eq. (I24l) . More realistically in larger networks a given isotope may be found in many reaction 
groups, some in equilibrium and some not, and for each group in equilibrium there will be a 
separate computed equilibrium value for the isotope. 

However, by hypothesis the equilibrium abundances of a given isotope computed for 
each of its equilibrated reaction groups cannot be too different, since this would violate the 
equilibrium assumption. For example, consider the a-network example from above. There 
is only one actual abundance Ya in the network at a given time, and its value must be such 
that it satisfies simultaneously the equilibrium conditions for reaction group I and reaction 
group II, within the tolerances of Eq. (l24l) : otherwise the equilibrium assumption would be 
invalidated. Therefore, restoration of equilibrium for a given isotope will correspond to setting 
its abundance to a compromise choice among each of the (similar) predicted equilibrium 
values for all equilibrated reaction groups in which it participates. This is a self-consistent 
approximation as long as the spread in possible equilibrium abundances remains consistent 
with the tolerances used to impose equilibrium in Eq. (I24l) . for if this were not true it would 
suggest that the equilibrium assumption itself is suspect. 

It is important conceptually that the use of "equilibrium assumption" as convenient 
shorthand in this discussion not be misinterpreted. The decision that a reaction group is in 
equilibrium is not an arbitrary choice but rather is made by the full network itself in each 
timestep. A reaction group is in equilibrium if it satisfies Eq. (|24l) for every species in the 
group, so the only arbitrary choices are the tolerances e,-. If the condition (l24l) is satisfied in a 
timestep for some reaction group, it is a statement by the network that in its current state the 
equilibration timescale for the reaction group has been found to be fast enough to maintain 
approximate equilibrium over the timestep. 
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6.2. Specific Methods for Restoring Equilibrium 

As noted above, in realistic networks we must have a systematic way to restore equilibrium 
at the end of numerical timesteps since those reaction in equilibrium at the beginning of the 
timestep will generally be disturbed slightly away from equilibrium during the timestep by 
those reactions that are not in equilibrium. This effect is tiny for one reaction timestep, but 
will accumulate to unacceptable error over many timesteps if we do not correct for it in each 
timestep. In this section we outline three potential approaches to this problem. Although these 
approaches are general, it probably is easiest to understand them in the context of a specific 
application so we shall illustrate with a 4-isotope alpha network having species a, ^^C, ^^O, 
and ^*^Ne, that we already illustrated schematically in 95.11 Using explicit notation for the 
alpha network we include the reactions 

0: 3a^ 1 : a + ^^C^^^O 

(47) 

2: ^i2(-_^^^ 20^6 3: a + i^O^^o^e, 

which we shall reference below in terms of the number on the left side for each reaction, and 
the differential equations (l34l)-(IT7]) governing the abundances then become 

Ya = -kfh^+li'^hu-ki'haYn + ki'h.e 

+ kfh^^ - k?-haY2o - kfhaYxe + ^^^20 (48) 
Yn = kfhl-k['^hn-kfhaYn+k['hi(,-kfh^2+k^^haY2Q (49) 
716 = k['haYn - fcf' - kfh^Yie + kfh2Q (50) 
Yio = kfhl2 - k?h^Y2Q + k^^haYie - -^{^^20, (51) 
where we use a notation Y{^^C) = Y12 and so on, for the abundances (defined in Eq. Q), 

in) (n) 

and the rate parameters k^ and k^ refer to forward and reverse rates, respectively, for the 
reactions pairs labeled by the numbers on the left sides in Eq. (|47l) . Thus kf'^ is the rate 
parameter for 3a — )• ^^C and ^r^^'* is the rate parameter for ^^Ne — )• a + ^^O. 

6.2.1. Using Iteration to Reimpose Equilibrium Abundance Ratios If we now impose 
equilibrium on the reaction a + ^^O ^ ^Ojsjg obtain the constraint 

ki^^ 

YaY\6 = "(3)" ^20 = KY20- (52) 

Inserting this constraint into Eqs. (I48l)-(|5T1) eliminates the last two terms in each of Eqs. (l48l) . 
(l50l) . and (|5TI) . giving the reduced equations 

Ya = -kf^h^+ki'^hn-k^/^YaYn + ki'hie+ki^h^2~ki^^YaY2o (53) 

712 = fcfVa - ki'^^Yn - k^f'^YaYn + ki'hie - /rf + fcJ'V«72o (54) 

Yi6 = ki^haYi2-ki^hie (55) 

Y20 = 4''^Y^2-kPYaY2o (56) 

This is the same number of equations as before, but these equations should now be less stiff 
than the original Eqs. (|48l)-(l5T1) because the terms eliminated by substituting the equilibrium 
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condition (|52l) are associated with the reaction pair that was first to come into equilibrium, 
precisely because it involves fast reactions. However, if we impose the equilibrium condition 
(|52|) at the beginning of a timestep, it generally will no longer be satisfied at the end of the 
timestep because the other reactions are not in equilibrium and will change the abundances of 
a, ^^O, and ^^Ne through non-equilibrium transitions. Thus, we must adjust the populations 
at the end of the timestep to reimpose the condition (|52|) . assuming the equilibrium condition 
to still be satisfied for that reaction pair. In addition to the algebraic constraint (|52l) . we have 
the condition 

Lx = = 4Ya + UYn + 16716 + 20720 = 1, (57) 

i 

imposed by the requirement that the sum of the mass fractions be unity (conservation of 
nucleon number; see Eq. Q). In contrast to the constraint ( |52l ). which is valid only if partial 
equilibrium conditions for the reaction a + '^O ^ ^^Ne are satisfied, the constraint (|57l) is 
a conservation law that must always be satisfied. Thus, to reimpose equilibrium we have 
available two conditions, Eqs. (l52l) and (|57l) . 

We assume equilibrium at the beginning of the timestep in the a + ^^O ^ ^*^Ne reaction 
(only); thus we advance the solution through a timestep At by solving the reduced equations 
(l53])-(l56l) using the explicit asymptotic algorithm [|T]|. Let us denote the (numerically 
computed) abundances at the end of this timestep by 7, = (7a,7i2,7i6,72o). At the beginning 
of the timestep the condition (|52|) is satisfied, by hypothesis, but at the end of the timestep 
in general 7a7i6 ^ KY2o- But the assumption that a + ^^O ^ ^*^Ne is in partial equilibrium 
implies (loosely) that the characteristic timescale T for this reaction group must be less than 
the integration timestep. Thus, we may assume that if we had taken short enough timesteps 
the reaction pair a + ^^O ^*^Ne would have brought the populations for 7^, 7i6, and 72o 
back into equilibrium by the end of the actual timestep taken, and we use the algebraic 
conditions given by Eqs. (|52l ) and (1571) to reimpose a + ^^O ^ ^^Ne equilibrium at the end of 
the timestep. Eqs. (l52l) and (l57l) allow us to write 

Fo = Ya~ ^720 = Fi=Yie- f 72o = 
F2 = Ya + 37i2 + 47i6 + 5720 - ^ = 

This may be written as the vector equation F = 0, which we can solve by Newton-Raphson 
iteration for the unknowns 7^, 7i6, and 72o, starting the iteration from computed values Y^, 
7i6, and 72o, and taking 7i2 to be fixed at the computed value. In terms of the Jacobian matrix 
J defined by 



dF 
dY 
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720 


/ 



(59) 



a Newton-Raphson iteration step then corresponds to choosing an initial vector Y , computing 
the corresponding values of F and J, solving the matrix equation 

J5Y = -F (60) 
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for the increment dY, and adding this increment to the original Y to get a corrected 
Y ^Y + 5Y. The corrected Y can then be used as the starting point for a second iteration, 
and so on until a convergence criterion is satisfied. 

Because it involves a matrix solution, this method has the potential to spoil the linear 
scaling with network size that is attractive about explicit methods unless the matrix equations 
can be solved by means other than brute force. Because the matrices for large networks are 
expected to be sparse and well-conditioned, it is likely that solution methods giving scaling 
not too different from linear are possible, but this needs to be demonstrated for this method. 
For our tests we shall solve these matrix equations using standard matrix packages. 



6.2.2. Use Iteration to Reimpose Equilibrium Abundances As an alternative to restoring 
abundance ratios, we may seek to reimpose equilibrium by requiring that the individual 
abundances of all isotopes participating in equilibrium be set to their equilibrium values (|2TI) 
at the end of their numerical timestep. Suppose that at some integration timestep there are 
n reaction groups in equilibrium. There will be some number of isotopes participating in 
these reaction groups, with some isotopes participating in more than one equilibrated reaction 
group. Define a vector Y of the distinct isotopes that are participating in at least one 
equilibrated reaction group. If we assume that both a + ^^C ^ ^^O and a + ^^O ^ ^*^Ne are 
in equilibrium for the network (|47] ). N = 4 and n = 2, and the components of Y are 

Y, = {YaJl2jl6j20} (61) 

where = y(^2c) and so on. Now define a vector F of normalized differences between the 
value of Yi at the end of the numerical timestep and its equilibrium value calculated from (f2T)) 
at the end of the numerical timestep, and an entry imposing conservation of particle number, 
and require that it vanish. For the above example, we obtain the vector equation F = 0, with 
the components of F given by 

p. _ ^12-4'^ ^le-y/e' ra-YH^ ^16-4'' ^20-4'^ v _ i \ r^?^ 

~ 1 ' ' ' y(3) ' y(3) ' y{i) ' V"^'' 

I. "cc "n '16 'a 'id "20 ) 

where 7, denotes the actual abundance of species i at the end of the numerical timestep 
and the computed equilibrium value of species / at the end of the timestep is Y^^\ with 



the index i labeling the position in the vector (1611) and the index j labeling the reaction 
group for which Eq. (f2T)) is solved. The first three entries in the vector (|62|) are associated 
with restoring equilibrium in reaction group 1 and the next three entries are associated with 
restoring equilibrium in reaction group 3. Notice that the abundances Ya and Y\(y appear more 
than once (twice each, for this example) in the vector F. This is because these isotopes appear 
in reaction group 1 and in reaction group 3, both of which are assumed to be in equilibrium. 
For the last entry Y.x — 1 is defined by Eq. (|57l) . and is a requirement that total nucleon number 
be conserved in the equilibrium restoration step. 

Thus restoration of equilibrium at the end of the timestep, subject to conservation of 
particle number, corresponds to solving F{Y) =0 for the vector Y (note that this vector 
contains only those isotopes that are part of reaction groups in equilibrium, not all of the 
isotopes in the network). As outlined in the previous section, the equation F(Y) = can be 
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solved iteratively for Y by choosing an initial guess for Y, computing F{Y), solving the 
matrix equation (|60l) for the increment dY, computing the improved Y ^Y + 5Y, and then 
repeating until a convergence tolerance is satisfied. Since the equilibrium abundances for the 
isotopes are likely to change very slowly in successive timesteps, this Jacobian matrix may 
be expected to be almost constant between two successive timesteps, and it will not change in 
successive Newton-Raphson iterations since the equilibrium abundances are constant within 
a given timestep. 

As for the method described in the previous section, a Newton-Raphson iteration on 
matrix equations is employed, so it will be highly desirable to solve these equations by means 
that retain the nearly linear scaling of the explicit method. From the discussion we see that for 
large networks with many reaction groups in equilibrium the matrices will be extremely sparse 
and only slightly changed in successive timesteps (and unchanged in successive Newton- 
Raphson iteration steps within a given timestep). We may expect that this structure lends 
itself to solutions with favorable scaling behavior in large networks, but that has not been 
investigated yet. 



6.2.3. Reimpose Equilibrium Abundances by Averaging The methods described in the 
previous two sections for restoring equilibrium at the end of numerical timesteps suffer from 
a certain level of redundancy because we have seen that in partial equilibrium the isotopic 
abundances in a reaction group are not independent but evolve according to a single timescale 
given by Eq. (l23l) (see the general discussion in ^4.11) . Thus, within a single reaction group 
specification of the equilibrium abundance of any one isotope, or of the progress variable 
associated with the reaction group, specifies the equilibrium abundance of all species in the 
group. Furthermore, within a single reaction group the evolution of the species in the group 
to equilibrium naturally conserves particle number, by virtue of constraints such as those of 



Eq. (1151) that are listed for all five reaction group classes in Appendix Appendix A 



Let us exploit this by working with the progress variable A,- from each reaction group. If 
all the reaction groups were independent, then we could restore equilibrium for the progress 
variables for n reaction groups in equilibrium simply by requiring 

/ Ai-Ai \ 
A2 — A2 







(63) 



where A, denotes the equilibrium value of A/ computed from Eq. (|2T)) and relations like 
Eq. (I22I) . Once the equilibrium value of A, (or the abundance of any one of the isotopes 
in the reaction group) is computed, the equilibrium values for all other isotopes in the group 
than follow from constraints like Eq. (|22|) that are tabulated for all reaction group classes in 
Appendix [Appendix A[ No probability conservation constraint is included in Eq. (|63l) because 



each reaction group considered in isolation conserves particle number automatically in the 
evolution to equilibrium. The solution of Eq. (1631) is then trivial, consisting of setting A, = A; 
for all /. 
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The simple considerations of the preceding paragraph cannot be used in the form 
presented: the reaction groups are generally not independent, because isotopes of the network 
appear typically in more than one reaction group, as we have discussed in ^6. 1.31 The 
equilibrium condition implies that the equilibrium abundances of a specific isotope computed 
for each of the reaction groups in which it is a member must be similar, but exact equality 
generally will not hold because of the finite tolerance e, used to test for equilibrium in Eq. (l24l) . 
Thus, we restore equilibrium for each isotope participating in partial equilibrium at the end 
of a timestep by replacing its computed abundance with its equilibrium value averaged over 
all equilibrated reaction groups in which it participates. If the species / is a member of M 
equilibrated reaction groups, its restored equilibrium value at the end of the timestep is 

j 

where the equilibrium abundances Y^^^ are computed for each of the reaction groups from 
Eq. (EB. 

There is one further issue: although the evolution to equilibrium for individual reaction 
groups conserves particle number, the averaging procedure of Eq. (|64|) over multiple reaction 
groups will introduce a fluctuation in the particle number since the chosen value of % 
computed from the average will generally differ from the individual P/"'^ that were computed 
conserving particle number. This difference will be very small in any one timestep, but can 
accumulate to unacceptable failure of particle number conservation for a long integration. 
Thus, for each timestep, after equilibrium has been restored through Eq. (|64l) we recompute 
the total particle number in the network (now summing over all isotopes, whether participating 
in equilibrium or not), compare it with the total particle number at the beginning of the 
timestep, and rescale all 7, in the network by a multiplicative factor that restores the particle 
number to its initial value. 

We have shown in a variety of comparisons that this procedure for restoring equilibrium 
works very well, giving results (isotopic abundances and timestepping profile) that are 
essentially identical to those of the iterative matrix algorithms discussed in the previous two 
sections. The advantage of the present approach is simplicity and automatically-favorable 
scaling with network size. The restoration of equilibrium involves only a summation over 
reaction groups to compute averages, followed by a multiplicative renormalization by a scalar. 
Thus there are no matrix inversions and no Newton-Raphson iterations, so this algorithm for 
restoring equilibrium is very simple to implement and should scale linearly with the size of 
the network, straight out of the box. All of the following examples were computed using this 
averaging algorithm. 

7. A Simple Adaptive Timestepper 

For testing the algorithms described here an adaptive timestepper has been employed that 
is described in more detail in Refs. [[B |2l [El. This timestepper is far from optimized but 
it leads to stable and accurate results for the varied astrophysical networks that have been 
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tested. Thus it is adequate for our primary task here, which is to establish whether explicit 
partial equilibrium methods can even compete with implicit methods for stiff networks. 

8. Toy Models for Partial Equilibrium 

As a first illustration of using the partial equilibrium methods developed in the preceding 
sections, we take a very simple thermonuclear network including only three isotopes, ^He 
(cc-particles), ^^C, and ^^O, connected by the reactions 3a ^ ^^C and a + ^^C ^ ^^O. Thus 
we have a total of four reactions in two reaction groups, with 3 a # ^^C belonging to reaction 
group class C and a + ^^C ^^O belonging to reaction group class B. The partial equilibrium 



formulas for both cases have been summarized in Appendix A This is an extremely simple 
network but we shall see that it demonstrates in a rather transparent way most of the essential 
features of a partial equilibrium calculation. 

Figure|4]illustrates a simulation with this network at a constant temperature of 5 x 10^ K 
and constant density 1 x 10^ g cm^^. This calculation employs the explicit asymptotic method 
n], but tries to impose partial equilibrium if the abundances of all reactions in a reaction group 
are within 0.01 of their equilibrium abundances. The dash-dotted blue curve labeled 1 /Rmax 
estimates the maximum stable fully-explicit timestep as the inverse of the current fastest rate 
in the network. This curve is set by the ^^C — )• 3a reaction initially, but by the a + ^^C — )• ^^O 
reaction after log? ~ — 7.5. 

First consider a purely asymptotic approximation. At the beginning the maximum 
accurate timestep is smaller than the maximum stable explicit timestep, so the algorithm 
takes explicit timesteps with a size set by accuracy and not stability considerations (which 
we arbitrarily cap for this example at dt ~ 0.1? to ensure accuracy). Near log? = —5.2 the 
explicit timestep becomes equal to the maximum stable explicit timestep (intersection of dash- 
dotted blue and dotted black curves), which is decreasing with time at this point because the 
rate for a + ^^C — t- ^^O is increasing with time. No isotopes are yet asymptotic, so the explicit 
timestep begins to decrease in order to remain below the dash-dotted blue curve and thus 
maintain stability, with the fluctuations in the timestep curve representing the attempts by the 
timestepping algorithm to increase the timestep being thwarted by the stiffness instability. But 
around log? = —3.5 one of the isotopes (^^C) becomes asymptotic and the timestep begins to 
increase rapidly over the explicit stability limit. However, the explicit asymptotic algorithm is 
unable to take large enough timesteps to make more isotopes asymptotic and the asymptotic 
timestep saturates at late times near J? ~ 1 x 10^^ s (dotted green curve). 

On the other hand, for the asymptotic plus PE calculation we proceed as for the purely 
asymptotic calculation until at log? ~ —4.4 the a + ^^C ^ ^^O reaction is determined by the 
network to satisfy the partial equilibrium condition and the net flux from this pair of reactions 
is removed from the numerical integration by the PE algorithm. Because of the removal of 
these fast components, the maximum stable timestep for the PE calculation (dashed red curve 
in Fig.Stb), corresponding to the inverse of the fastest timescale remaining in the numerical 
integration) begins to increase relative to that for the asymptotic calculation (dash-dotted blue 
curve). In response, the PE integration step size (solid red curve in Fig. ID^b)) also begins to 
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Figure 4. The effect of imposing partial equilibrium on integration timesteps for the network 
a ^ '^C '^O at a constant temperature of 5 x lO'^ K and constant density 1 x 10^ g cm^^, 
assuming initial mass fractions of 0.50 for '^C and ^^O. Reaction rates were taken from 
standard compilations lfT2l . (a) Mass fractions as a function of time. The results are 
almost indistinguishable between the purely asymptotic calculation (Asy) and the asymptotic 
plus partial equilibrium calculation (PE). (b) The purely-explicit integration timestep (black 
dotted), the asymptotic timestep (green dotted), the partial equihbrium timestep (solid red), 
the approximate maximum explicit timestep for the asymptotic calculation (dash-dotted blue), 
and the maximum stable timestep for the partial equilibrium calculation (dashed red), (c) The 
fraction of isotopes that are asymptotic and the fraction of reaction groups that are equilibrated 
as a function of time. The dotted green curves in the middle and bottom figures refer to the 
results one obtains with the explicit asymptotic method if partial equilibrium is not imposed. 
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Log time (s) Log time (s) 

Figure 5. (a) Equilibration and reaction timescales for Fig. |4] (b) Stiffness stability domains 
implied by the equilibration timescales in part (a). 

increase until at log? ~ —3.4 it is about 100 times larger than the purely asymptotic timestep. 
At this point, the network determines that the 3a ^^C reaction group also satisfies the partial 
equilibrium condition and removes this pair of reactions from the numerical integration too. 
Now, since there are only two equilibrium reaction groups in our simple network and they 
have both been removed by the PE algorithm, the numerical integrator is effectively solving 
the set of equations dY = 0. Since all timescales have now been removed from the numerical 
integration by the PE algorithm, there is no stiffness instability at all and the maximum stable 
timestep (dashed red curve) goes to infinity. 

Provided that the two reaction groups remain in equilibrium, we are now free to take 
timesteps limited only by accuracy. Accordingly the PE timestep increases rapidly and once 
again becomes equal to an upper limit of dt ~ O.lt imposed by our arbitrary overall accuracy 
constraint (solid red curve). Thus, after a network evolution time of one second the partial 
equilibrium timestep is approximately four orders of magnitude larger than that for the explicit 
asymptotic algorithm, with the calculated isotopic abundances as a function of time being 
indistinguishable in the two cases. This translates into a total integration time that is 400 times 
faster for the asymptotic plus partial equilibrium calculation versus the purely asymptotic 
calculation in the example of Fig. HI 

A deeper understanding of the factors that determine the timesteps in Fig. |4] may be 
obtained by consulting Fig. [51 which displays all timescales relevant for the problem as a 
function of time. Figure [Sj^a) shows the timescales defined by the inverse of the reaction rates 
for all reactions in the network (dotted and dashed curves), the timescales T of Eq. (|23l) for 
establishing equilibrium for the two reaction groups in this calculation (magenta and blue 
bands), and the integration timestepping for various methods (solid curves). The differential 
equations to be solved are those of Eqs. (|48l) -(l5TI) with terms involving ^*^Ne removed: 

Ya = ~kf'h^+ki'^hu-4'^YaYn + ki'h,e (65) 
712 = kf^h^-ki'^hn-ki'haYn+ki'hie (66) 
Yi6 = ki'haYn - ki'hie - kfhaY.e (67) 
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Figure 6. Equilibration of the single reaction group 3a ^ C. (a) Mass fractions in 
asymptotic and asymptotic plus partial equilibrium approximation, (b) Timescales and 
integration timesteps. The dashed curves are the timescales set by the inverses of the reaction 
rates; the dotted red curve is the equilibration timescale (|23] ) for the reaction group. The solid 
green curve is the integration timestepping for an asymptotic approximation and the solid red 
curve is the timestepping for an asymptotic plus partial equiUbrium approximation. 



In the explicit asymptotic calculation (where we ignore the role of the equilibrium timescales 
denoted by the magenta and blue bands) we see that the timestep is rigidly limited by the 
inverse rate for a + ^^C — )■ ^^O (dashed purple curve) until ^^C becomes asymptotic (labeled 
^^C asy). This partially lifts the timestepping restriction because Eq. (|65l) and its source terms 
for a + are then removed from the numerical integration. However, since no 

other isotopes become asymptotic, the numerical integration is unable to get past the ceiling 
imposed by the timescales associated with ^^C — )■ 3a and ^^O — t- a + ^^C that remain in 
Eqs. (I66l)-(l67l) and the asymptotic timestep saturates around 10^ s at late times. 

In contrast, the partial equilibrium method is able to remove all timescales in all 
three equations associated with the reaction group a + ^^C ^ ^^O when the timestep dt is 
comparable to the magenta band denoting x{a + ^^C # ^^O). This complete removal of a set 
of fast timescales replaces the original problem with one that is less stiff. That permits the PE 
method to increase its timestep sufficiently that dt quickly reaches the equilibrium timescale 
associated with 3a ^ ^^C (blue band) and removes completely all timescales associated with 
this reaction group too, leaving an integration problem having no stiffness limitation on the 
explicit timestep. 

An even simpler picture of the relationship between asymptotic and partial equilibrium 
approximations emerges if we consider the evolution of a single reaction group. In Fig.[6ta), 
by evolving a single reaction group consisting of the triple-alpha reaction and its inverse, we 
see clearly the role of the individual reaction timescales and the timescale for approach to 
equilibrium of 3a ^^C in setting the timestepping for asymptotic integration. For accuracy 
in this example, we have arbitrarily limited the integration timestep to dt < 0. It. As illustrated 
in Fig.|6tb), initially, the largest stable fully-explicit timestep is governed by the inverse of the 
rate for the 3a — )■ ^^C reaction. Until logf ~ 4, the largest accurate timestep lies below this 
timescale, so there is no stiffness limitation on the explicit integration. At around log? ~ 4 
the ^He (a particle) becomes asymptotic and the asymptotic timestep begins to exceed the 
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limit set by the 3 a ^^C reaction by a small amount. However, the other isotope in the 
network, ^^C, never satisfies the asymptotic condition over the entire range of integration. 
Thus, the asymptotic method is never able to remove completely the stiffness timescales set 
by the 3a ^^C reactions (since they will remain at all times in the differential equation for 
^^C, which must be integrated numerically). Thus the asymptotic-method integration timestep 
saturates at late times near the timescale set by the inverse of the rate for 3a — t- ^^C. 

On the other hand, in the asymptotic plus PE calculation the reaction group 3a ^^C 
becomes equilibrated according to the criteria of Eq. (I24l) at about the time the equilibrium 
timescale for the reaction group (shown as the magenta band) crosses the dt curve near 
log? ~ 3. Since the reaction group is now assumed to be in equilibrium, all flux terms 
associated with both 3a — )■ ^^C and ^^C — )■ 3a are removed from the numerical integration and 
the numerical integrator must solve the system dYa/dt = dY^/dt = 0. Since all timescales 
have been removed from the network, the corresponding integration timestep is set only by 
accuracy criteria. As a result, the partial equilibrium integration in Fig.[6]is more than 300,000 
times faster than the asymptotic integration of the same system, with essentially identical 
results. 

Similar considerations apply in more complex networks and it is useful to view the 
timescales of Fig.^a) or Fig.[6tb) as establishing stiffness stability domains in the dt versus t 
plane that are illustrated in Fig.[5tb). Thus, a fully-explicit integration is restricted by stability 
considerations to the lower domain marked I, an asymptotic (or QSS) calculation is restricted 
to domains I and II, an asymptotic or QSS plus partial equilibrium calculation is restricted to 
domains I, II, and III, and all of these methods would be unstable in domain IV. However, as 
illustrated in Fig. Ob), if we are sufficiently clever we can push domain IV far enough into 
the upper left corner of the dt — t plane that it plays little practical role because the instability 
domain corresponds to timesteps that would be undesirable from an accuracy point of view. 

The examples displayed in Figs. [5] and [6] are not very complex and it should be not at all 
surprising that the description of the system becomes simple and the stable timesteps large 
at later times, since the mass fraction curves become almost constant in this region. But 
that is the whole point of the partial equilibrium approach: near equilibrium the system is 
complicated when viewed in terms of competing independent reactions, but becomes simple 
when viewed in terms of groups of reactions coming into equilibrium and thus bringing the 
whole system into equilibrium. These examples illustrate for some very simple networks, but 
the same principles will apply to the more complex networks to which we now turn. 

9. Tests on Some Thermonuclear Alpha Networks 

We have tested the partial equilibrium algorithm described in earlier sections in a variety of 
thermonuclear alpha networks. In this section we give some representative examples of those 
calculations. Because they are extremely challenging reaction network problems to solve, 
we shall concentrate on examples representative of conditions expected in Type la supernova 
simulations (temperatures in the range 10^-10^" K, densities in the range 10^-10^ gcm^-^, 
and initial equal mass fractions of ^^C and ^^O). Such conditions lead quickly to high degrees 
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Table 4. Reactions of the alpha network for partial equilibrium calculations. The reverse 
reactions such as ^"Ne a + '^O are photodisintegration reactions, 7+ ^"Ne a + ^^O, 
with the photon 7 suppressed in the notation. 
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of equilibration in the approach to QSE (quasi-statistical equilibrium) and NSE (nuclear 
statistical equilibrium), and are a stringent test of partial equilibration methods. Our long- 
term goal is application of these methods to larger networks, but an alpha network provides a 
highly-stiff test system that is small enough to allow significant insight into how the algorithm 
functions. We also note as a practical matter that the most ambitious published calculations 
for thermonuclear networks coupled to hydrodynamical simulations have employed alpha 
networks. 

The reactions and corresponding reaction groups used in the calculations are displayed 
in Table IH We use the standard REACLIB library [fT2]| for all rates except for three of 
the heavy-ion reactions (corresponding to reaction groups 3, 5, and 6 in the table) which 
have inverse reactions that are absent from the REACLIB compilation and were taken 
from Ref. [fT4||. For typical Type la supernova conditions the rates for the capture and 
photodisintegration reactions in Table Hlbecome comparable and as large as ~ 10^^ — 10^^ s^^. 
The corresponding maximum stable fully-explicit integration timestep will be of order the 
inverse of the maximum rate in the network, so timesteps as short as ~ 10^^^ s may be 
required for stability of a fully-explicit integration. On the other hand, the characteristic 
physical timescale for the primary Type la explosion mechanism is of order one second. Thus 
integration of the alpha network coupled to a hydrodynamics simulation of a Type la explosion 
could require 10^^ or more explicit network integration timesteps, which is obviously not 
practical and indicates graphically that this is an extremely stiff system. 
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Figure 7. Partial equilibrium calculation for constant Tg = 5 and p = 1 x 10^ g cm^^ alpha 
network, (a) Mass fractions. Solid lines are implicit calculations (Xnet) and dashed lines 
are asymptotic + PE calculations, (b) Integration timesteps for several different integration 
methods. Magenta dash-dot is the implicit code Xnet ifTSll . black dashed is the semi-implicit 
code YASS |16| (reproduced from from Ref. ifTOl ). solid red is the current asymptotic 
plus partial equilibrium result, solid black is the QSS plus partial equilibrium calculation 
reproduced from Ref. 1 10|, black dash-dot is a QSS calculation using the formalism of Ref. 
12], dotted green is an asymptotic calculation using the formalism in Ref. 1 1 1, and dashed blue 
is the estimated maximum stable fully-explicit timestep. 



9.1. Comparisons of Explicit and Implicit Integration Speeds 

We shall be comparing explicit and implicit methods using codes that are at very different 
stages of development and optimization. Thus, they cannot simply be compared directly 
with each other. We assume that for codes at similar levels of optimization the primary 
difference between explicit and implicit methods would be in the extra time spent in implicit- 
method matrix operations. Hence, if the fraction of time spent on linear algebra is / for 
an implicit code, we assume that an explicit code at a similar level of optimization could 
compute a timestep faster by a factor of F = 1 / ( 1 — /) . Factors F have been tabulated in 
Ref. [IJ for several networks. Then we may compare roughly the speed of explicit versus 
implicit codes (possibly at different levels of optimization) by multiplying F by the ratio of 
implicit to explicit integration steps required for a given problem. This procedure has obvious 
uncertainties, and likely underestimates the speed of an optimized explicit versus optimized 
implicit code, but will give a useful lower limit on how fast the explicit calculation can be 
relative to implicit methods. 

9.2. Constant Intermediate-Temperature, Low-Density Example 

A calculation for constant Tg = 5 and density of 1 x 10^ g cm^-^ in an alpha network is shown 
in Fig.|71 This network has 16 isotopes with 48 reactions connecting them, and a total of 19 
reaction groups that can separately come into equilibrium. For this and all other calculations 
equilibrium conditions were imposed using Eq. (l24l) . with a constant value e,- = 0.01 for 
the tolerance parameter. The calculated asymptotic plus partial equilibrium (Asy + PE) 
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mass fractions are compared with those of a standard implicit code in Fig. [Tfa). There are 
small discrepancies in localized regions, especially for some of the weaker populations, but 
for the most part the agreement is rather good. Although we shall generally display mass 
fractions down to 10^^^ in our examples for reference purposes, it is important to note that 
for reaction networks coupled to hydrodynamics only larger mass fractions (those larger than 
say ~ 10^^ — 10^^) are likely to have significant influence on the hydrodynamics. Thus, 
the largest discrepancies between mass fractions calculated by PE and implicit methods in 
Fig. Ha) mostly imply uncertainties in the total mass being evolved by the network of less 
than one part in a million, which would be completely irrelevant in a coupled hydrodynamical 
simulation. 

Timestepping for various integration methods is illustrated in Fig. Hb), where there are 
several interesting features to discuss. 

(i) The PE timestepping is seen to compare rather favorably with standard fully-implicit 
and semi-implicit codes. The PE calculation required 3941 integration steps while 
the implicit code Xnet required only 600 steps, but this factor of 6.5 timestepping 
advantage is offset substantially by the expectation that for a 16-isotope network an 
explicit calculation should be ~ 3 times faster than the implicit calculation for each 
timestep [1]. Thus we conclude that for fully-optimized codes the implicit calculation 
would be perhaps several times faster for this example. The semi-implicit timestepping 
curve was reproduced from another reference but a comparison of the curves in Fig. |7tb) 
suggests that an optimized partial equilibrium code is likely to be at least as fast as the 
semi-implicit YASS code for this example. 

(ii) The Asy -i- PE and implicit timestepping are seen to be many orders of magnitude better 
at late times than the results of our purely asymptotic (Asy) and purely quasi-steady- 
state (QSS) calculations, which are based on the formalisms discussed in Refs. [HI 13. 
The reason is apparent from Fig. [8l For times later than about log? = —6, significant 
numbers of reaction groups come into partial equilibrium, which asymptotic and QSS 
methods are not designed to handle. 

(iii) The earlier application by Mott et al fTO^ of a QSS plus partial equilibrium calculation 
for this system lags orders of magnitude behind the current implementation of asymptotic 
plus partial equilibrium in timestepping at late times. In fact, the timestepping from Ref. 
[[Toll is only a little better than that of the pure QSS results from the present paper. 

In Fig. [9] we display the partial equilibrium timescales of Eq. (l23l) that are associated with the 
calculation of Fig.|71 with the reactions that are removed by partial equilibrium or asymptotic 
considerations indicated (compare with the simpler earlier example in Fig. [S]). In Fig. [9ta) 
each reaction has been removed from the graph when it reaches equilibration. In Fig. |9tb) 
we also indicate in dashed gray those reactions that have been at least partially removed from 
the numerical integration by isotopes becoming asymptotic, or ranges of times for reactions 
in which reactants have very small concentration and therefore make only small contribution. 
The complete removal of reactions by partial equilibration and the partial removal of reactions 
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Figure 8. Fraction of reactions that are treated as being in equilibrium as a function of 
integration time for the calculation of Fig. [T] (solid red curve). Also shown is the fraction of 
isotopes that become asymptotic in the PE calculation (dotted green curve). 




Figure 9. Timestepping and partial equilibrium timescales associated with the calculation 
in Fig. [T] (a) The equilibrium timescales, with the reactions removed for times after they 
become equilibrated in the PE calculation, (b) As for part (a), but with labels suppressed 
and with time ranges for which equilibrium reactions have been partially removed by the 
asymptotic approximation indicated by dashed curves. The heavy, dark curves indicate times 
for which the corresponding equilibrium timescale has been removed neither by asymptotic 
nor PE approximations and therefore is fully operative in producing stiffness. These curves 
clearly have a large influence on the maximum timestep dtp^ that the PE calculation is able to 
take. 




Figure 10. Partial equilibrium calculation for constant Tg =7 and p = 1 x 10^ g cm^^ alpha 
network, (a) Mass fractions. Solid curves are asymptotic + PE, dashed curves are asymptotic, 
and dotted curves are implicit (Xnet). (b) Integration timesteps for the present asymptotic plus 
partial equilibrium calculation (solid red), the implicit code Xnet [15] (dotted blue), and the 
asymptotic method of Ref. UJ (dashed green). 



by asymptotic conditions clearly play a significant role in determining the envelope of the 
maximum stable and accurate partial equilibrium timestep, as seen in Fig. Hb). 

9.3. Constant Higher-Temperature, Intermediate-Density Example 

A calculation for constant Tg = 7 and density of 1 x 10^ g cm^^ in an alpha network is shown 
in Fig. \Wi Again we see extremely good agreement among the mass fractions calculated 
by an implicit method, the asymptotic method, and the asymptotic plus partial equilibrium 
method. The timesteps for the partial equilibrium and implicit methods are quite similar, with 
the implicit method requiring 510 integration steps and the partial equilibrium method 506 
integrations steps. An optimized explicit code can probably evaluate a timestep ~ 3 times 
faster than an implicit code for a 16-isotope network (T\ and we conclude that for similarly- 
optimized codes the partial equilibrium calculation would be slightly faster. In contrast, the 
asymptotic method lags far behind at late times, because the equilibration becomes significant 
after 10^^ s, with the fraction of reactions that are equilibrated reaching 100% for times later 
than 10^4 s. As a consequence, the integration of Fig. \T0\ required more than two million 
purely-asymptotic integration steps. 

9.4. Example with a Hydrodynamical Profile 

The preceding two examples employed alpha networks at extreme but constant temperature 
and density. A partial equilibrium calculation using a hydrodynamical profile with the 
dramatic temperature rise characterizing a thermonuclear runaway in a Type la supernova 
simulation is illustrated in Fig. [TTJ The thermonuclear runaway in degenerate white-dwarf 
matter subjects the reaction network to extreme conditions in this example: the difference 
between the fastest and slowest reaction timescales is more than 1 1 orders of magnitude and 
during the thermonuclear flash the temperature increases by 6.6 billion K in only 1.5 x 10^ 
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Figure 11. Asymptotic, PE, and implicit calculations with an alpha network for Type la 
supernova conditions, (a) Mass fractions, (b) Integration timesteps. (c) Hydrodynamic profile 
ifTTl . (d) Fraction of isotopes asymptotic and partiaUy-equilibrated reactions. 



s (a rate of temperature change corresponding to 4.4 x 10 K/s). The asymptotic plus PE 
method required 712 steps to integrate this problem, compared with 524 steps for the implicit 
calculation, suggesting that an optimized explicit code would be about twice as fast as the 
implicit code because it can compute a timestep about three times faster. In contrast, the 
purely asymptotic calculation gave accurate results but required almost 190,000 integration 
steps for the time range shown. 

As may be seen from Fig.fTlTd). the partial equilibrium calculation becomes completely 
equilibrated after log? ~ —5.2 and after that time the PE timesteps are limited only by 
accuracy. Choosing to restrict the subsequent timesteps to dt = 0.3t, continuation of the 
PE (or implicit) integration of Fig. [TT] until the characteristic timescale for the supernova 
explosion of ~ 1 s requires only a few tens of additional integration steps. In contrast, if the 
asymptotic integrator continued to take the timesteps it is taking at the end of the calculation 
in Fig. [TTJ it would require more than a billion additional integration steps to reach the full 
physical timescale for the supernova. 

9.5. Synopsis 

The three examples shown above are representative of a number of problems that we 
have investigated with the new asymptotic plus partial equilibrium methods. Our general 
conclusion is that for alpha networks under the extreme conditions corresponding to a Type 
la supernova explosion the asymptotic plus partial equilibrium algorithm is capable of giving 
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timestepping that is orders of magnitude better than asymptotic or QSS approximations alone 
when the system approaches equilibrium, that these timesteps are often competitive with 
current implicit and semi-implicit codes, and that they are orders of magnitude faster than 
previous attempts to apply such explicit methods to extremely stiff thermonuclear networks. 

10. Improving the Speed of Explicit Codes 

The preceding discussion, and that of Refs. |2l, have assumed that we can compare 
(approximately) the relative speed of explicit and implicit methods that are presently 
implemented with different levels of optimization by comparing the number of integration 
timesteps required to do a particular problem with each method. The assumption of this 
comparison is that once the matrix overhead associated with implicit methods is factored out, 
the rest of the work required in a timestep (computing rates and fluxes, and so on) is similar 
for explicit and implicit methods. While this is useful for a first estimation, it should not 
be viewed yet as a quantitative comparison of the speed of implicit and explicit methods for 
two basic reasons. (1) The current timestep routines for implicit methods and the explicit 
methods presented here are at very different levels of algorithmic optimization. (2) Because 
the timestep routines are not yet equivalently optimized, we have not made a systematic study 
of optimal timestepping parameters for either the implicit or explicit methods. Thus, the 
results presented here should be interpreted as indicating that for the problems investigated 
explicit methods can take stable and accurate timesteps that are comparable to those of implicit 
methods. Further work will be required to determine whether implicit or explicit methods are 
faster and by how much for specific problems. 

There are at least three reasons why the integration speed for explicit methods may 
potentially be even faster than that estimated by such a procedure in this paper and in 
Refs. Eia. 

First, in this paper we have emphasized partial equilibrium methods in league with 
asymptotic approximations. But in Ref. [2] we presented evidence that quasi-steady-state 
(QSS) methods give results comparable to asymptotic methods, but often with timesteps that 
can be as much as an order of magnitude larger. Thus, we may expect that a QSS plus partial 
equilibrium approximation could give even better timestepping than the examples shown here. 
We have not yet investigated this possibility systematically. 

Second, we have seen that the timestepping algorithm employed in this paper and in 
Refs. [|Il|2l is not very optimized. On the one hand, results such as those exhibited in Fig.|9tb) 
suggest that our timestepper is tracing qualitatively the explicit-integration stability boundary. 
But that same figure suggests that in many regions explicit timesteps as much as an order of 
magnitude larger could still be stable. Thus, it is possible that we have over-estimated the 
number of integration steps required by our explicit methods in the examples discussed in this 
paper because of an adequate but not optimized timestepper. This possibility also remains to 
be investigated. 

Third, and potentially of most importance, our demonstration that algebraically- 
stabilized explicit methods may in fact be viable for large, extremely-stiff networks implies 
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new considerations with respect optimization. For standard implicit methods, the most 
effective optimizations involve improving the numerical linear algebra, since in large 
networks more than 90% of the computing time for implicit methods can be consumed 
by matrix inversions. But with algebraically-stabilized explicit methods the bulk of the 
computing time for large networks is spent in computing the rates. Thus, the most effective 
optimization for our explicit methods lies in increasing the speed with which reaction rates 
are computed. For example, can effective schemes to compute rates that exploit modem 
multi-core architectures with GPU accelerators be developed for realistic simulations of 
reaction networks coupled to hydrodynamics? Such optimizations could be used for implicit 
integrations too, but they will be less effective in increasing their speeds because in large 
networks determining reaction rates is only a small part of the computing budget of an 
implicit method. In contrast, we would expect our explicit integration methods to gain speed 
substantially with faster methods to compute reaction rates. 

11. Extension to Larger Thermonuclear Networks 

As we have shown elsewhere, well away from equilibrium asymptotic and quasi- steady- state 
approximations can compete with or even out-perform current implicit codes in a variety of 
extremely stiff networks 111121. There are various important problems, such as the nova and 
tidal supernova examples discussed in Refs. lUlHl, where the system does not equilibrate 
strongly. Even in problems where equilibrium is important overall, there will be many 
hydrodynamical zones and timesteps for which equilibrium is not very important. For these 
cases, the asymptotic and QSS approximations alone may be viable. However, to compete 
across the board with implicit solvers, the present paper makes clear that asymptotic and QSS 
approximations must be supplemented by partial equilibrium methods to remain viable in 
those zones where non-trivial equilibrium processes are present. 

The examples of the previous section have shown that asymptotic plus partial equilibrium 
methods can solve stiff thermonuclear alpha networks with timestepping in the same ballpark 
as implicit or semi-implicit codes, and accuracy sufficient for coupling to fluid dynamics 
codes. Although alpha networks represent the current state of the art in coupling reaction 
networks to hydrodynamics in such problems, our primary goal is to use these methods 
to extend network sizes to physically more reasonable values (hundreds or thousands of 
isotopes). Because the explicit methods that we are developing can compute a timestep 
faster, and scale approximately linearly and therefore more favorable with network size than 
implicit methods, their relative advantage grows with the size of the network. Thus, for 
explicit methods to realize their full potential we must extend the previous partial equilibrium 
examples to include very stiff networks containing hundreds of reacting species. That work 
is underway and will be reported in future publications, but in this section we make some 
general remarks and report some preliminary results concerning these efforts. 
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Figure 12. (a) Maximum rates as a function of time in an alpha network and a 150-isotope 
network for constant density and a hydrodynamic temperature profile given in part (b). Rates 
were computed using the compilation in Ref. [12|. (b) Hydrodynamical temperature profile 
corresponding to the calculation shown in part (a); the density is approximately constant over 
this time range. This profile is characteristic of evolution under Type la supernova conditions 
in a single zone of a carbon-oxygen white dwarf with an initial = 3 and p = 1 x 10^ g cm^^. 
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11.1. Relative Stiffness of Alpha and Realistic Thermonuclear Networks 

The primary difference (other than size) between an alpha network and a more realistic 
thermonuclear network is that a more realistic network will typically contain protons and 
neutrons. The rates for reactions with protons and neutrons are often many orders of 
magnitude higher than those induced by alpha particles; as a result, realistic networks tend 
to contain much broader ranges of characteristic timescales than alpha networks, and thus 
are often much stiffer, as illustrated in Fig. [I2l The fastest rates in the realistic network are 
seen to be 6-8 orders of magnitude larger than those in the alpha network for this example. 
These very fast reactions will tend to bring reaction groups into equilibrium rapidly. Thus, 
we find that realistic networks are much stiffer than alpha networks, but the most important 
effect for an explicit integration is not the increased magnitude of the stiffness itself (which 
the asymptotic and QSS approximations can handle in a manner competitive with implicit 
methods, as documented in Refs. [[Ill2l), but rather that these fast reactions quickly manifest 
themselves in fast equilibration timescales, which asymptotic and QSS algorithms cannot deal 
effectively with unless supplemented by a partial equilibrium approximation. 

In Fig. \T3\ we compare how quickly reactions come into equilibrium (a) for an alpha 
network and for a 150-isotope network under constant temperature and density conditions 
typical of some zones in a Type la supernova explosion, and (b) using a hydrodynamical 
profile derived from one zone of a Type la supernova simulation. In this calculation we 
have not yet implemented the partial equilibrium approximation, but have integrated using 
the asymptotic approximation and tracked the number of reactions that would satisfy the 
partial equilibrium condition as a function of time in that case. We see that the effect is quite 
pronounced, with the large, realistic network having a much greater fraction of its reactions 
satisfying the partial equilibrium condition from very early times when compared with the 




Figure 13. Fraction of equilibrated reactions as a function of time for an alpha network (dotted 
blue, labeled a) and a 150-isotope network (solid red, labeled 150). (a) Constant temperature 
Tg = 5 and constant density p = 10^ gcm^^ case, (b) Case with hydrodynamic profile having 
density essentially constant and the temperature variation given by the dashed green line, with 
Tg indicated on the right axis. In both cases we see that the 150-isotope network equilibrates 
much faster and much more extensively than the alpha network. In these examples we have 
assumed a reaction group to be in equilibrium if all isotopes participating in the reactions of 
the group have abundances y, within 1% of their equilibrium values. 



alpha network. For example, by the time the alpha network begins to show any substantial 
equilibration the 150-isotope network is 30-40% equilibrated in both cases. This suggests 
that the partial equilibrium approximation may have an even more dramatic effect on the 
speed with which large networks can be integrated than we have found for alpha networks. 

11.2. A Toy Model: Protons Plus Alpha-Particles 

Let us illustrate the ideas of the previous section with a toy model that is simple to understand, 
but that contains many of the essential features expected for a realistic network. In Fig. [14] we 
display a calculation with a 4-isotope alpha network containing the isotopes ^He, ^^C, ^^O, 
and ^*^Ne, including all reactions among these species from the REACLIB [12] library. We see 
that the partial equilibrium abundances are essentially identical to those calculated using only 
the explicit asymptotic algorithm, and that at the end of the calculation the partial equilibrium 
calculation is taking timesteps that are about 3 orders of magnitude larger than those for the 
explicit asymptotic method alone, and about 4 orders of magnitude larger than the maximum 
timestep that would be stable for a standard explicit method. 

In Fig. \T5\ we repeat the calculation of Fig. [I4l except that we add to the 4-isotope 
alpha network two new isotopes: protons and ^^F. Now the network contains new reactions 
involving protons that are much faster than the alpha reactions and hence is much stiffer. 
Again we find that the mass fractions calculated using the asymptotic plus partial equilibrium 
method and the asymptotic method alone are almost identical, but the disparity in the 
corresponding integration speeds is now much larger. At the end of the calculation illustrated 
in Fig. [151 we see that the PE integrator is taking timesteps that are about 10^ times larger 
than those for the explicit asymptotic integrator, and almost 10^ times larger than would be 
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Figure 14. 4-isotope alpha network for constant temperature Tg = 5 and density p = 
1 X 10^ gcm^^. Reactions rates were taken from Ref. lfT2]| and initially there were equal 
mass fractions of '^C and ^''O. (a) Mass fractions (asymptotic solid; asymptotic plus partial 
equilibrium dotted), (b) Timestepping (solid red for asymptotic plus partial equilibrium and 
dotted green for asymptotic). The maximum stable purely-explicit timestep is estimated by 
1 //?max and shown as the dashed blue curve. 
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Figure 15. 4-isotope alpha network plus protons and F, for constant temperature = 5 
and density p = 1 x 10** gcm^^^. Reactions rates were taken from Ref. [12| and the initial 
mass fractions were Xq(}^C, '*'0, 'H, ^^F) = (0.500, 0.400, 0.001,0.122). (a) Mass fractions 
(asymptotic solid; asymptotic plus partial equilibrium dotted), (b) Timestepping (solid red 
for asymptotic plus partial equilibrium and dotted green for asymptotic). The maximum 
stable purely-explicit timestep is estimated by 1/^max and shown as the dashed blue Une. 
It is determined by the rate for '^F p + ^^O across the entire time range shown. This 
photodisintegration rate is constant because the temperature is constant. 



Stable with a normal explicit method without the partial equilibrium approximation. 

This simple example illustrates succinctly the arguments of the previous section. 
Comparing Fig. [15] with Fig. O we see that the addition of a single set of proton reactions 
has an enormous influence on the stiffness of the network. In the pure a-particle case of 
Fig. [14l the maximum stable fully-explicit timestep in the approach to equilibrium is of 
order 10^^ s and is set by a-particle reactions. In this case a sufficient number of isotopes 
become asymptotic to permit an asymptotic timestep of near 10^ seconds in the approach to 
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equilibrium. This is not a competitive timestep, since we see from Fig.[14]that the PE method 
is able to take timesteps about three orders of magnitude larger, but it is far better than the 
situation in Fig. \T5\ There we see that the maximum stable fully-explicit timestep over the 
entire range of integration is of order only 10^^-^ s, because it is set by proton reactions that 
are much faster than any a-particle reaction, as we have already illustrated in Fig.[l2l (Under 
these conditions the fastest of these is the photodisintegration 7+ ^^F — t- p + ^^O.) 

In Fig. \T5\ the PE method is taking timesteps in the approach to equilibrium that are 
six orders of magnitude larger than those of the pure asymptotic method. The reason for 
this increased disparity between the PE plus asymptotic and pure asymptotic methods is not 
primarily a better timestep for the PE method (it is taking timesteps in both calculations that 
are near the limit imposed by accuracy constraints), but rather the large depression of the 
maximum stable asymptotic timestep produced by adding the much faster proton reactions 
to the network. However, we see that when the fast proton reactions are removed from the 
numerical integration by the PE approximation, the timestep of Fig. [15] becomes comparable 
to that for the much less stiff pure alpha network of Fig. [H 

Hence, our toy model supports rather well the speculations of the previous section. 
Realistic large networks are much stiffer than a-networks, but this increased stiffness is 
associated largely with fast neutron and proton reactions that may be more susceptible to 
improvement by partial equilibrium approximations than simple alpha networks. The results 
of this section suggest that in large realistic networks the systematic removal of the fastest 
timescales by the PE approximation may permit timestepping by these explicit methods that 
is comparable to that of implicit methods, even when equilibrium is important. As we have 
seen, timestepping comparable to that of implicit methods in large networks implies that an 
explicit partial equilibrium code should be considerably faster than the implicit code, because 
it can compute each timestep faster. 

12. Summary and Conclusions 

Numerical integration for extremely stiff reaction networks has typically been viewed as 
requiring implicit schemes, because normally explicit timesteps of efficient size would be 
unstable. An alternative to implicit integration modifies the equations using approximate 
solutions to reduce the stiffness, and then integrates the resulting equations explicitly. For 
example, asymptotic and quasi-steady-state methods found a measure of success integrating 
moderately stiff chemical-kinetics systems explicitly EdHKlOl, but have generally failed for 
extremely stiff systems, giving inaccurate results with non-competitive integration timesteps. 
In this paper, and the two preceding it ID 13, we have presented evidence suggesting that 
algebraically- stabilized explicit integration can give correct results and timesteps competitive 
with that of implicit methods for even the stiffest of networks. 

Central to this new view of the efficacy of explicit integration for stiff systems is the 
understanding that in reaction networks there are at least three different sources of stiffness: 

(i) Negative populations, which can occur for an initially small positive population if an 
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explicit timestep is too large; this anomalous negative population leads to exponentially 
growing terms that destabilize the network. 

(ii) Macroscopic equilibration, where the right sides of the differential equations dY = F = 

— F approach a constant derived from the difference of two large numbers (the 
total flux in and total flux out F ), and numerical errors in taking this difference 
destabilize the network if the timestep is too large. 

(iii) Microscopic equilibration, where on the right sides of Eq. ([T]) the net flux in specific 
forward-reverse reaction pairs (/-^ — /j-^) tends to zero as the system approaches 
equilibrium; this leads to large errors if the timestep is too large because the net flux 
is derived from the difference of two large numbers and the timescale equilibrating the 
populations is short compared with the desired numerical timestep. 

The algebra required to stabilize the explicit solution depends on which of these forms of 
stiffness is operative. In Refs. [[Tl 111 we have shown that asymptotic and quasi-steady- 
state approximations give correct results, with timestepping highly competitive with implicit 
methods even for the stiffest of thermonuclear networks, as long as the system exhibits only 
the first two types of stiffness listed above. However, in those papers we have also shown 
that asymptotic and quasi-steady-state methods give correct results, but with timesteps that 
are much larger than for purely explicit methods but still far from competitive with implicit 
methods, if there is any significant microscopic equilibration in the system. 

In this paper we have addressed in depth the role of microscopic equilibration in stiffness 
and presented a partial equilibrium scheme that can be used to augment asymptotic and quasi- 
steady- state methods when equilibrium becomes important in a network. We have presented 
strong evidence that this partial equilibrium scheme can remove the timestepping deficiencies 
encountered by asymptotic and quasi-steady-state methods in the approach to equilibrium, 
making these methods competitive with implicit methods even near equilibrium. The methods 
that we have presented build on earlier work |l3l [HI [TOj, but we find that our versions of 
asymptotic and quasi-steady-state methods are far more accurate, and our versions of partial 
equilibrium are faster by multiple orders of magnitude, than found in previous work when 
applied to extremely stiff astrophysical thermonuclear networks. 

The present paper, in conjunction with the previous ones on asymptotic methods 
and quasi- steady- state methods 0, suggests that algebraically-stabilized explicit integration 
methods are capable of stable timesteps competitive with those of implicit methods even 
for extremely- stiff reaction networks. Explicit methods can compute timesteps faster 
than implicit methods in large networks, implying that algebraically- stabilized explicit 
algorithms may compete favorably with implicit integration, even in complex, extremely- 
stiff applications. Because explicit methods scale linearly with network size, these findings 
could permit the coupling of physically more realistic reaction networks to fluid dynamics 
simulations in a variety of fields. 
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Appendix A. Appendix: Reaction Group Classification 



Applying the principles discussed ^4.51 to the reaction group classes in Table [2] gives the 
following systematic partial equilibrium properties of reaction group classes for astrophysical 
thermonuclear networks. 

Reaction Group Class A (a ^ b) 

Source term: ^ = + k,yb Constraints : ya+yb = ci = yl-^yl 
Equation: ^=bya + c b = -kf c = k. Solution: ya{t) = y^^' - f (l - 
Equil. solution: = ~f = Equil. timescale: = ^ = ^ 

Equil. tests: ^^y^ < (z = a,b) Equil. constraint: = 
Other variables: = c\—ya 

Progress variable: X=yl-ya = - A Jz, = + A 

Reaction Group Class B (a + b ^ c) 

Source term: ^ = -hyayb + Kyc 

Constraints: = ci = y\-y^a yb +yc = C2=yh + yc 

Equation: ^ = ayl + bya+c a = —kf b = —{cikf + kh) c = k^{c2 — Ci) 

Solution: Eq. (fT9l ) Equil. solution: Eq. (|2TI) Equil. timescale: Eq. (|23] ) 

Equil. tests: ^-^j^ < (z = a,b, c) Equil. constraint: = ^ 

Other variables: = ci +3^0 yc = C2—yb 

Progress variable: X=y'^-ya = yl] - A = - A = 3;° + A 
Reaction Group Class C (a + b + c ^ d) 

Source term: ^ = -kfyaytyc + Kyd Constraints: ya - = ci = j° - 3;° 

I (ya +yb+yc)+yd = c3 = ^{y^,+yt+y^c)+ yd y^ - yc = c2 = yl-y'^ 

Equation: ^ =ayl + by a + c 

a = -%^ + fcf(ci +C2) b= -{kiClC2+h) C = (C3 + ^Cl + ^C2)fcr 

Solution: Eq. (fT9l) Equil. solution: Eq. (|2TI) Equil. timescale: Eq. (|23l ) 
Equil. tests: ^^^^ < £/ (/ = a,b,c,d) Equil. constraint: ™^ = | 

Other variables: = - ci yc=ya-C2 yd = C3-ya + ^{ci+ C2) 
Progress variable: A = j° - j« y^=y'^-X yi,=y^-X y,.=y^.-X 

yd = ^+yd 
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Reaction Group Class D (a + b ^ c + d) 

Source term: ^ = -kfjayb + hycya Constraints: - y^, = ci = y° - y^ 

ya+yc = c2 = yl+y^c ya+yd = ci=yl+y^d 

Equation: = ay^ + by a + c a = k^ — ki b = —kr{c2 + C3) + k^ci c = kj^C2C3 

Solution: Eq. (fT9l) Equil. solution: Eq. (|2TI) Equil. timescale: Eq. (|23l) 

Equil. tests: ^-^^ < £, (/ = a,b,c,d) Equil. constraint: JgJ = | 

Other variables: = - ci j,. = C2 - ya yd = 03- ya 

Progress variable: X=y^^-ya = - A 3;^ = - A = + A 

Reaction Group Class E(a + bv^c + d + e) 

Source term: ^ = -kfyayb + krycyaye 

Constraints: ya + \{yc+yd +ye) =ci=yl + \{y^, +3^0 

Equation: ^ = ayl + by a + c 

a = (3ci -j^)fcr-^f ^ = C2fcf- (a/3 + a7+j87)A:r c = ha^y 

a = Ci + ^(c3+C4) /3 = Ci - |C3 + ^C4 7=Ci + ^C3-|c4 

Solution: Eq. (fT9l ) Equil. solution: Eq. (|2TI) Equil. timescale: Eq. (|23]) 
Equil. tests: ^^C^ < £, (/ = a, b, c, J, e) Equil. constraint: ^ggj = | 
Other variables: yb=ya-C2 yc = oc-ya yd=p-ya ye = y-ya 
Progress variable: A = j° - j« 3;^ = - A 3;^ = 3;° - A 3;^ = 3;° + A 
Jrf = Jrf + A = y° + A 



In the equilibrium tests we have allowed the possibility of a different e, for each species /, but 
in practice one would often choose the same small value e for all The results presented in 
this paper typically have used = 0.01 for all species. 
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